Computer simulation of model cohesive powders: 
influence of assembling procedure and contact laws on low consolidation states. 
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Molecular dynamics simulations are used to investigate the structure and mechanical properties of 
a simple two-dimensional model of a cohesive granular material. Intergranular forces involve elastic- 
ity, Coulomb friction and a short range attraction akin to the van der Waals force in powders. The 
effects of rolling resistance (RR) at intergranular contacts are also studied. The microstructure of 
the cohesive packing under low pressure is shown to depend sensitively on the assembling procedure 
which is applied to the initially isolated particles of a granular gas. While a direct compression 
produces a final equilibrated configuration with a similar density to that of cohesionless systems, 
the formation of large aggregates prior to the application of an external pressure results in much 
looser stable packings. A crucial state variable is the ratio P* — Pa/Fo of apphed pressure P, 
acting on grains of diameter a, to maximum tensile contact force Fq. At low P* the force-carrying 
structure and force distribution are sensitive to the level of velocity fiuctuations in the early stages 
of cluster aggregation. The coordination number of packings with RR approaches 2 in the limit of 
low initial velocities or large rolling friction. In general the force network is composed of hyper- 
static clusters, typically comprising 4 to a few tens of grains, in which forces reach values of the 
order of Fq, joined by barely rigid arms, where contact forces are very small. Under growing P*, it 
quickly rearranges into force chain-like patterns that are more familiar in dense systems. Density 
correlations are interpreted in terms of a fractal structure, up to a characteristic correlation length 
^ of the order of ten particle diameters for the studied solid fractions. The fractal dimension in 
systems with RR coincides, within measurement uncertainties, with the ballistic aggregation result, 
in spite of a possibly different connectivity, but is apparently higher without RR. Possible effects of 
micromechanical and assembling process parameters on mechanical strength of packings are evoked. 

PACS numbers: 

81.05.Rm: Porous materials, granular materials, 

83.10.Rs: Computer simulation of molecular and particle dynamics, 
61.43.Hv: Fractals, macroscopic aggregates, 
47.57.J-: Colloidal systems 



I. INTRODUCTION 

Granular materials are currently being studied by 
many research groups [E 0, H, motivated by funda- 
mental issues (such as the relations between microstruc- 
ture and global properties) as well as by practical needs 
in civil engineering and in the food and drug industries. 
The relation of their mechanical behavior in quasistatic 
conditions to the packing geometry, which depends itself 
on the assembling procedure, tends to escape intuition 
and familiar modelling schemes. 

The configuration of the contact networks is hardly 
accessible to experiments, even though particle positions 
are sometimes measured [^, 0, @|, and some experi- 
mental quantitative studies on intergranular contacts car- 
ried out in favorable cases (such as millimeter-sized beads 
joined by capillary menisci 0, 0, [13, [HI ) • Intergranular 
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forces are also, most often, inaccessible to measurements. 
Consequently, computer simulation methods of the "dis- 
crete element" type, as introduced 30 years ago have 
proved a valuable tool to investigate the internal states 
of granular systems. Simulation methods like molecular 
dynamics [Tj] or "contact dynamics" [l^ have been 
gaining an increasingly large constituency of users and a 
wide range of applications, as witnessed e.g., by recent 
conference proceedings 

Dry assemblies of grains interacting via contact elas- 
ticity and friction, such as sands or glass beads, might 
form stable packings of varying solid fraction (typically 
between 58% and 64% for monosizcd spheres if they do 
not crystallize), which deform plastically in response to 
changes in stress direction, rather than stress intensity. 
Their elastic or elastoplastic properties have been stud- 
ied by discrete simulation (see, e.g., [l^ . [l3|)j and, in 
agreement with laboratory experiments and macroscopic 
modelling [l^ , found to depend sensitively on the initial 
density. Numerical simulation also stressed the impor- 
tance of additional variables such as coordination num- 
ber [l9| and fabric [13, HH , and it has often been applied 
to the study of quasi-static stress-strain behavior of gran- 
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ular assemblies (refs. [T^ [13, are a few examples 

among a large literature). 

Cohesive grains exhibit much larger variations in their 
equilibrium densities, and they are sensitive to stress in- 
tensity as well as direction : on increasing the confining 
pressure, the specific volume of a clay can irreversibly 
decrease by a factor of 4 [23|. Likewise, series of ex- 



periments carried out in the Seville group on model pow- 
ders [iSl (xerographic materials) in which the strength of 
van der Waals attraction is controlled by additives cov- 
ering part of the grain surfaces, reveal a similar variation 
of porosity with confining pressure. It is notable that 
such packings of particles of rotund shape and nearly the 
same size can stay in mechanical equilibrium at much 
lower solid fractions (down to 25-30%) than cohesionless 
granular systems. 

Despite this wider variety of equilibrium structures and 
mechanical behaviors, cohesive granular materials have 
much less frequently been investigated by numerical sim- 
ulation than cohesionless ones. 

Some of the recent numerical studies, such as those 
of refs. [2^, [2^ have investigated the packing struc- 
tures of spherical beads deposited under gravity, depend- 
ing on micromechanical parameters, including adhesion 
strength. Another set of publications report on simula- 
tions of the dynamical collapse and compaction, both in 
two [27I . [2M , [2g | and three [sO] dimensions, the main re- 
sults being the relations between density and pressure in- 
crements, and their dependence on micromechanical pa- 
rameters. Some works focussed on the fracture of bound 
particle assemblies in static [sil [s^ or dynamic [s^ con- 
ditions, others on wet bead packs in which cohesion stems 
from liquid bridges joining neighboring particles, investi- 
gating the structure of poured samples [sj] or the shear 
strength [Tll| of such materials. These two latter types 
of studies deal with relatively dense materials, as does 
the numerical biaxial compression test of [s^. Flow of 
cohesive materials has also been addressed in recent pub- 
lications 36, 37, 38]. 

Yet, numerical studies of the mechanics of loose, solid- 
like cohesive granulates are quite scarce. This contrasts 
with the abundant literature on the geometry of model 
loose particle packings and colloidal aggregates, which 
tend to form fractal structures. Refs. [3^, |40| are use- 
ful overviews of aggregation processes and the geometric 
properties of the resulting clusters, as obtained by nu- 
merical simulation. In such processes, particle aggregates 
are usually regarded as irreversibly bound, rigid solids, 
while the interaction between separate clusters reduces to 
a "sticking rule" , so that both intra- and inter-aggregate 
mechanical modelling is bypassed. Interestingly, one sim- 
ulation study j4l| shows that structures resulting from 
geometric deposition algorithms are not always stable 
once a mechanical model is introduced. 

It seems that numerical simulations of both geometric 
and mechanical properties of loose granular assemblies 
forming solid aggregates are still lacking. 

The present paper addresses part of this issue. It re- 



ports on numerical simulation studies of cohesive granu- 
lar materials, with the following specificities: 

• the assembling process is simulated with the same 
mechanical model as applied to solid-like configura- 
tions, and its influence on the packing microstruc- 
ture is assessed ; 

• special attention is paid to loose particle packings 
in eguilibrium under vanishing or low applied pres- 
sure; 

• both geometric and mechanical properties are in- 
vestigated ; 

• isotropic and homogeneous systems are studied, as 
representative samples for bulk material properties. 

We consider a simple model system in two dimensions, 
introduced in section [TTl along with the numerical simu- 
lation procedure. Despite its simplicity we shall see that 
this model yields results that are amenable to compar- 
isons with experimental situations. 

Section. IIIII is devoted to the important issue of the 
procedure to prepare samples, and its influence, as well 
as that of micromechanical features such as rolling re- 
sistance (RR), on final density and coordination num- 
ber in solid packings in equilibrium. In Section IIVI we 
investigate the force distributions and force patterns of 
the equilibrated loose configurations under vanishing or 
low applied pressure. Some specific aspects of the force- 
carrying structures in low density assemblies will be stud- 
ied and related to the assembhng process. In Section IVl 
we characterize the geometry and density correlations in 
loose samples, resorting to the fractal model traditionally 
employed for colloidal aggregates. Finally we conclude 
in Section IVII with a few remarks about future improve- 
ments and further developments of this work, some of 
which will be presented in a forthcoming publication [i^] . 



II. MODEL MATERIAL 
A. System definition, equations of motion 

We consider a two-dimensional model material: an as- 
sembly of N disks with diameters {di)i<i<N uniformly 
distributed between a/2 and a. The maximum diameter, 
a, will be used as unit of length. The mass of grain i is 
rui = df/a'^ and its moment of inertia li = midf/8, i.e. 
disks are regarded as homogeneous bodies and the mass 
of a disk of maximum diameter a is the unit of mass. 

The disks are enclosed in a rectangular cell the edges 
of which are parallel to the axes of coordinates xi and 
X2, with respective lengths Li and L2. Periodic bound- 
ary conditions are used, thereby avoiding wall effects. 
Neighboring grains, say i and j, might interact if they are 
brought into contact or very close to each other, hence 
a force Fij and a moment Fy exerted by i onto j at the 
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contact point. Simulations do not model material defor- 
mation in a contact region, but consider overlapping par- 
ticles, and the contact point is defined as the center of the 
intersecting surface of the two disks. In the case of an in- 
teraction without contact, the force will be normal to the 
surfaces at the points of nearest approach, and therefore 
carried by the line of centers. Let denote the position 
of the center of disk i. rij = rj — ri is the vector joining 



the centers of i and j , and hi 



{di + dj)l2 their 



overlap distance. The degrees of freedom, in addition to 
the positions r^, are the angles of rotation 9i, velocities 
Vi, angular velocities LOi — 6i of the grains (1 < « < N), 
the dimensions {La)a=i.2 of the cell containing the grains 
and their time derivatives, through the strain rates: 

in which denotes the initial size for the corresponding 
compression process. The time evolution of those degrees 
of freedom is governed by the following equations. 
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In Eqns. ^ and only those disks j interacting with i, 
i. e. in contact or very close, will contribute to the sums 
on the right-hand side. In Eqn. ([3|), is the externally 
imposed stress component, cr*^ is the measured stress 
component, resulting from ballistic momentum transport 
and from the set of intergranular forces Fij, A — L1L2 
denotes the cell surface area, and M is a generalized in- 
ertia parameter. 

Stresses an and a22, rather than strains or cell dimen- 
sions, are controlled in our simulation procedure. Note 
that compressions are counted positively for both stresses 
and strains. Eqn. ([3]) entails that the sample will expand 
(respectively, shrink) along direction a if the correspond- 
ing stress cr^ is larger (resp., smaller) than the requested 
value (T^jj , which should be reached once the system equi- 
librates. This barostatic method is adapted from the ones 
initially proposed by Parrinello and Rahman [H, [3, 
for Hamiltonian, molecular systems. 

The choice of the "generalized mass" M is rather arbi- 
trary, yet innocuous provided calculations are restricted 
to small strain rates. In practice we strive to approach 
mechanical equilibrium states with good accuracy, and 
choose M in order to achieve this goal within affordable 



computation times. We usually attribute to M a value 
equal to a fraction of the sum of grain masses (3/10 in 
most calculations), divided by a linear size L of the cell. 
This choice is dimensionally correct and corresponds to 
the appropriate time scale for strain fluctuations in the 
case of a thermodynamic system. 



B. Interaction law 

The contact law in a granular material is the relation- 
ship between the relative motion of two contacting grains 
and the contact force. As we deal with particles that may 
attract one another at short distance without touching, 
the law governing intergranular forces and moments is 
best referred to simply as the interaction law. 

Although the interaction we adopted is based on the 
classical linear "spring-dashpot" model with Coulomb 
friction for contact elasticity, viscous dissipation and slid- 
ing, as used in many discrete simulations of granular me- 
dia [13, [13, [H, \id\ , some of its features (short-range at- 
traction and rolling resistance) are less common ; more- 
over, one can think of different implementations of the 
Coulomb condition, depending on which parts of normal 
and tangential force components are taken into account. 
Therefore, for the sake of clarity and completeness, we 
give a full, self-contained presentation of the interaction 
law below. 

We express intergranular forces in a mobile system 
of coordinates with axes oriented along the normal unit 
vector hij (along r^ ) and the tangential unit vector iy 
(3) {^ijT^ij is ^ direct base in the plane), and use the con- 
vention that repulsive forces are positive. 

The intergranular force Fij , exerted by grain i onto its 
neighbor is split into its normal and tangential compo- 
Tijtij thus defining scalars Nij and 



nents, Fij — Nijfiij 
Tij. Nij comprises a static term depending on the dis- 
tance between disk centers, combining contact elasticity 
and distant, van der Waals type attraction, as shown on 
Fig. [T^, and a velocity-dependent viscous term N^^ . Tij 
(Fig. [T|d) is due to the tangential elasticity in the con- 
tact, and is limited by the Coulomb condition. If disks i 
and j are not in contact, both the tangential component 
of force Fij and the viscous part of the normal compo- 
nent vanish, while i and j still attract each other if the 
gap {hij > 0) between their surfaces is smaller than the 
attraction range Dq (0 < hij < Dq): 



with 



^Fo (1 



-) n^j (4) 



This expression is a linear approximation of a realistic 
van der Waals force law (see Fig.[ljL), and contains two es- 
sential parameters, maximum attractive force Fq, range, 
Dq. Typically, Fq is of the order of 7?, 7 b eing a superfi- 
cial energy, I the typical size of asperities [43| and Do is 
in the nanometer range. 

In the case of contacting disks {hij < 0), the attractive 
term N"', is kept constant, equal to —Fq, while strains in 
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FIG. 1: Graphical representation of the model for the adhesive elastic contact force as a function of the distance between 
the surfaces of particles i and j, hij. (a) The elastic normal force consists of a repulsive Hookean part N^^ plus a linearized 
attractive part N°'j. (b) The elastic tangential force is limited by the Coulomb cone (adhesion shifting its tip to — -Fb on the 
normal force axis). 



the contact region result in normal (Nf^) and tangential The elastic tangential force in contact i,j is linearly 
{Tij) elastic forces. It is also assumed that a viscous 



related to the elastic part Suf^ of the total relative tan- 



normal term N^j opposes relative normal displacements. 
One thus writes: 



- Tij tij 



(5) 



The different terms introduced in Eqn. (O are defined 
according to the following models. First, 



Nf: 



-KmH 



Nllij 



is the linear elastic unilateral repulsion, due to the normal 
deflection —hij in the contact as the disks are pressed 
against each other. is the normal stiffness coefficient, 
related to the elastic moduli of the material the grains 
are made of. 

The viscous normal force opposes the normal relative 
receding velocity Svfj = hij ■ [vj — Vi) as long as the 
contact persists. The relative normal motion of two disks 
i and j in contact is that of an oscillator with viscous 
damping, and rjij is the damping coefficient. We choose 
its value as a constant fraction ( of the critical damping 
coefficient. 



I AK^minij 



(6) 



This is equivalent to the choice of a constant restitution 
coefficient in normal collisions if fo = 0. In the presence 
of attractive forces the apparent restitution coefficient in 
a collision will depend on the initial relative velocity, and 
will be equal to zero for small values, when the receding 
velocity after the collision will not be able to overcome 
the attraction and separate the particles. The minimum 
receding velocity for two particles of unit mass [i.e., of 
maximum diameter a) to separate is V*^/2, with 



V* = ^FoDo. 



(7) 



gential displacement Auf, , as 



and is subject to the Coulomb inequality. Kt is the 
tangential stiffness coefficient. Au^- can be updated for 
all closed contacts according to 



dt 



[Vij ■ tij) 



and vanishes as soon as the contact opens. Its elastic 
part satisfies 



dSuf 
dt 



= H 



I^Nfj 
Kt 



\5uiA (% -ty) 



in which H denotes the Heaviside function. This last 
equation introduces the friction coefficient /i. It is im- 
portant to note that the Coulomb inequality. 



(8) 



applies to the sole repulsive elastic component of the nor- 
mal force (see Fig. [TJ)). We chose not to implement any 
tangential viscous force. 

The moment that disk i exerts onto its contacting 
neighbor j, of radius Rj, in its center, is denoted by Fy 
in Eqn. It is first due to the tangential contact force, 
then to a possible moment F^J of the force density dis- 
tribution within the contact region. One thus writes: 



(9) 



FjJ is most often neglected on dealing with smooth, con- 
vex particle shapes, because the contact region is very 
small on the scale of the particle radius. 
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To model rolling resistance (RR) , like in [48*1 , we intro- 
duce a rotational stiffness parameter and a rotational 
friction parameter in contacts, so that rolling elas- 
ticity and rolling friction are modelled just like sliding 
elasticity and friction. One thus writes 



r'- 



Kr- 59,. 



while enforcing the inequality 

Kr\5e,,\<^irN^J. 

This involves the definition of 5t 
of the total relative rotation AOij. The total relative 
rotation angle satisfies 



(10) 

as the elastic part 



dt ^ 
while the equation for 59ij is 



dt 



H 



Kr 



|(5% I {ujj - Ui). 



Parameters and Hr are often related to the size of a 
contact region [23| ■ K^. is dimensionally the product of a 
stiffness by the square of a length, which is of the order of 
the contact size. In the following we set Kr to 10~'^a?KM, 
while /ir, which has the dimension of a length, is chosen 
equal to 10~^/ia. The motivation for the introduction of 
RR into our model is twofold. First, cohesive particles 
are usually small (typically less than 30/xm in size) and 
irregular in shape. Contacts between grains are likely 
to involve several asperities, and hence some lateral ex- 
tension, of the order of the distance between asperities, 
however small the normal deflection —h. Then, it will 
be observed that even quite a small rotational friction 
has a notable influence on the microstructure of cohesive 
packings. 



C. Control parameters and dimensional analysis 

In this section we present the dimensionless param- 
eters which express the relative importance of different 
physical phenomena. Such parameters enable qualitative 
comparisons with real materials, bearing in mind that the 
present model is admittedly an idealization of real pow- 
ders and that our simulations do not aim at quantitative 
accuracy. 

Dimensionless numbers related to contact behavior are 
the reduced interaction range Do/a, the friction coeffi- 
cient fi, the viscous damping parameter C, and the stiff- 
ness parameter k. 

Under the attractive force —Fq, the elastic deflection 
of one contact is 



ho = Fo/K, 



N 



(11) 



The stiffness parameter n = aKj^/Fo characterizes the 
amount of elastic deflection ho under contact force Fo, 



relative to grain size a {ho /a = k A suitable analo- 
gous deflnition for Hertzian spheres in three dimensions 
would he K = [Ea'^/Fof/^. 

The dimensionless number ho /Do is the ratio of 
elastic to adhesive stiffnesses, and its physical mean- 
ing is similar to that of the Tabor parameter A — 
{l/Do){'j'^a/E^y^^ [49j] for a Hertzian contact between 
spheres of diameter a when the material Young modu- 
lus is E and the interfacial energy is 7 (more precisely, 
the equilibrium normal deflection ho, due to adhesion, in 
the contact between an isolated pair of grains, satisfies 
A ^ [ho/DoY^^ in this case). 

The viscous damping parameter, C,, corresponds to a 
normal restitution coefficient e^ = exp[— tt^/ ^Jl — C/^] in 
the absence of cohesion {Fo = 0). 

In our calculations we set C, — 0.8, corresponding to 
a high viscous dissipation in collisions, or a very low 
restitution coefficient in binary collisions. Models with 
a constant C were adopted in other published simulation 
works m, [s^l , although little is known about dissipation 
in collisions. C is known to influence the packingstruc- 
tures obtained in the initial assembling stage [4^ [5l| . 
but we did not investigate its effects in the present study. 
The simulations reported i n |25| . [2^ use the viscous force 
model introduced in ref. |52l ]. with a choice of param- 
eters corresponding to strongly overdamped dynamics 
{i.e., analogous to ^ 1 in our case). 

In addition to those control parameters determined by 
the contact behavior, other dimensionless numbers are 
introduced by the loading or the process being applied 
to the material. The effect of the external pressure, com- 
pared to the adhesion strength, is characterized by a di- 
mensionless reduced pressure P* : 



P* = Pa/Fo 



(12) 



In the present paper, we focus on the assembling process 
and the low P* range. As we shall see below (Sec. IIIip 
low density, tenuous structures are then stabilized by ad- 
hesion, and the relevant force scale is Fo. However, as 
briefly reported in js^l, such structures tend to collapse 
upon increasing P* . These phenomena will be the sub- 
ject of another paper Wolf et al. [1^ introduced a 
dimensionless stress proportional to P* , and observed, in 
numerical simulations, stepwise increases in pressure to 
produce large dynamical collapse effects around P* = 1. 
The importance of P* was also stressed in simulations of 
cohesive granular flow, in which the effects of cohesion 
on rheological laws were expressed in terms of a cohesion 
number defined as 1/P* [37|. In three dimensions, P* 
should be defined as a^P/Fo. 

For large reduced pressures, externally imposed forces 
dominate the adhesion strength, and one should observe 
behaviors similar to those of confined cohesionless gran- 
ular materials. For P* > 1, the relevant force scale is 
aP. The influence of k, which should then be defined as 
K = /P, so that the typical contact deflection h satis- 
fies h/a <x. \yas studied in simulations of grains with- 
out adhesion [5J|. Whatever the reference force used to 
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define it, the limit of rigid grains is k ^ +00. With rela- 
tively soft grains (say, k below 10'^), a significant number 
of additional contacts appear in dense configurations, due 
to the closing of gaps between near neighbors. Such a k 
parameter defined with reference to pressure, in the case 
of contacts ruled by Hertz's law between spherical grains 
made of a material with Young modulus E, should be 
chosen as k = (E'/P)^/'^, in order to maintain h/a ~ k,~^. 

In order to stay within the limit of rigid grains both 
for small and large P*, we choose quite a large value of 
K = KMa/Fo: k = 10'^ or k = 10^. 

Table U summarizes the values (or the range of val- 
ues) of dimensionless parameters in the simulations pre- 
sented below. In addition to those values of the param- 





c 


K 


Kt 


Do 


Kr 




P* 




a 


Kno^ 


a 


0.15, 0.5 


0.8 


lo^ lo"* 


1 




10-* 


0, 10-^ 


0, 0.01 



TABLE I: Values of dimensionless model parameters used in 

most simulations. Note that ho /Do is fixed by k = — and 

Fo 

Do/a to 10^^ or 10"^. In the absence of cohesion, or for 
values of P > Fo/a, K is defined as Km/P 



eters, adopted as a plausible choice for realistic orders 
of magnitudes, some calculations were also performed 
with deliberately extreme choices, such as very large RR 
(fir — 0.5a) or absence of friction (fi — and /i^ = 0), in 
order to better explore some connections between mi- 
cromechanics and macroscopic properties. The corre- 
sponding results will be described in Section HVl 

The definition of dimensionless parameters, suitably 
generalized to three-dimensional situations as P* = 
o^P/Fq and k ~ {Ea^ / FqY/^ for spherical particles of di- 
ameter a, enables one to discuss qualitative features and 
orders of magnitude in the model system defined with the 
parameters of Table |T] with comparisons to some cohesive 
packings studied in the laboratory. 

When adhesive forces are due to liquid menisci joining 
neighboring particles, we should take Pq ~ 7^*; where 7 is 
the surface tension. P* = 1 corresponds then to confining 
pressure P in the range of 10-100 Pa for millimeter-sized 
particles, taking standard values for 7. Those are rather 
low pressures in practice, which are comparable, e.g., to 
the ones caused by the weight of a typical laboratory 
sand sample. Thus wet granular materials are commonly 
under reduced pressures P* of order 1 or larger, and are 
not observed with much lower solid fractions than dry 
ones @,i,|iS[il|. 

The cohesive powders studied in refs. [l^, [5^ [5^ [s^l 
are xerographic toners with typical particle diameter 
a ~ 10 ^m. Pq, the van der Waals attractive force, is 
a few tens of nN, and the range Dq is several nanome- 
ters [11]. Therefore, a reduced pressure P* = 0.01 would 
correspond to about 1 Pa in the experimental situation 
[sgj . This is an initial state of very low consolidation 
stress, which is present in a powder under gravity, pro- 



vided a controlled gas flow, going upwards through the 
powder, counterbalances part of its weight [59J. As to 
contact stiffnesses, our values of h^/a would correspond 
to S - 0.1 GPa (for Kn = 10'*Po/a) or 3.2 GPa (for 
Kn — lO^Po/a); while the ratio Dg/a would imply an in- 
teraction range of 10 nm. This gives us the correct orders 
of magnitudes for the toner particles, those being made of 
a relatively soft solid (polymer, such as polystyrene) with 
E^3—6 GPa. Xerographic toner particle s appea r to un- 
dergo plastic deformation in the contacts j56l. [60l. [6ll. [6^ . 
Plastic deflections of contacts are accounted for in the 
model of ref . (ssj , applied to the simulation of a biaxial 
compression of a dense powder. In our study, for simplic- 
ity's sake, and because we expect macroscopic plasticity 
of loose samples to be essentially related to the collapse 
of tenuous structures, we ignored this feature. 

D. Equilibrated states 

Although numerical simulations of the quasistatic 
response of granular materials requires by definition 
that configurations of mechanical equilibrium should be 
reached, equilibrium criteria are sometimes left unspec- 
ified, or quite vaguely stated in the literature. Yet, in 
order to report results on important, often studied quan- 
tities like the coordination number or the force distribu- 
tion, it is essential to know which pairs of grains are in 
contact and which are not. Due to the frequent occur- 
rence of small contact force values, this requires forces 
to balance with sufficient accuracy. We found that the 
following criteria allowed us to identify the force-carrying 
structure clearly enough. We use the typical intergran- 
ular force value Pi = max(Po,Pa) to set the tolerance 
levels. A configuration is deemed equilibrated when the 
following conditions are fulfilled: 

• the net force on each disk is less than 10^'' Pi, and 
the total moment is lower than 10^"* Pi a; 

• the difference between imposed and measured pres- 
sure is less than 10"^Pi/a; 

• the kinetic energy per grain is less than 5 • lO^^Pia. 

We observed that once samples were equilibrated accord- 
ing to those criteria, then the Coulomb criterion ([5]), as 
well as the rolling friction condition (jTU]) were satisfied as 
strict inequalities in all contacts. No contact is ready to 
yield in sliding, and with RR no contact is ready to yield 
in rolling either. 

III. ASSEMBLING PROCEDURE 

It has been noted in experiments [2^ and simula- 
tions [H, [H, that the internal structure and resulting 
behavior of solid-like granular materials is sensitive to the 
sample preparation procedure, even in the cohesionless 
case. 
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In the case of powders, it has been observed that the 
sedimentation in dry nitrogen (to minimize the capiUary 
effects of the humidity on the interparticle adhesion) of 
a previously fluidized bed produced reproducible states 
of low solid fractions (down to 10 — 15%) (63i. i64|. This 
initial state under such a low consolidation, as we com- 
mented in III C[ plays a decisive role on the evolution of 
the dynamics of powder packing. That is, appreciable 
differences in initial states will lead to considerable ones 
in final packings W^. This is mainly due to the role of 
aggregation, which we shall analyze in the second part of 
this section. 

The motivation of this section is to investigate the de- 
pendence on packing procedure in a cohesive granular 
system, the first step being to obtain stable equilibrated 
configurations with low densities. For comparison, some 
simulation results are presented for the same model ma- 
terial with no cohesion. 

Specimens were prepared in two different ways, respec- 
tively denoted as method 1 and method 2, and the result- 
ing states are classified as type 1 or type 2 configurations 
accordingly. 

Due to our choice of boundary conditions, our sam- 
ples will be completely homogeneous, under a uniform 
(isotropic) state of stress. This choice is justified by 
the complexity of seemingly more "realistic" processes, 
such as gravity deposition, due to the influence of many 
material (such as viscous dissipation, as recalled in Sec- 
tion lll (J|l and process parameters. Both pouring rate and 
height of free fall should be kept constant during such 
a pluviation process in order to obtain a homogeneous 
packing [m], [6^1 with cohesionless grains. Cohesive ones, 
because of the irreversible compaction they undergo on 
increasing the pressure, would end up with a density in- 
creasing with depth. Hence our choice to ignore gravity 
in our simulations. Our final configurations should be 
regarded as representative of the local state of a larger 
system, corresponding to a local value of the confining 
stress. 



A. Method 1 

In simulations of cohesionless granular materials, a 
common procedure [H, [13, [111 to prepare solid samples 
consists in compressing an initially loose configuration (a 
"granular gas"), without intergranular contacts, until a 
state of mechanical equilibrium is reached in which in- 
terparticle forces balance the external pressure (further 
compaction being prevented by the jamming of the parti- 
cle assembly). We first adopted this traditional method, 
hereafter referred to as method 1, to assemble cohesive 
particles. 

In this procedure, disks are initially placed in random 
non-overlapping positions in the cell, with zero velocity. 
We denote such an initial situation as the /-state. Then 
the external pressure is applied, causing the cell to shrink 
homogeneously. Thus contacts gradually appear and the 



configuration rearranges until the system equilibrates at 
a higher density. 

Examples of equilibrated configurations are shown on 
Fig. [21 with and without cohesion. This state is char- 
acterized by its solid fraction (<& — J^i'^'^l f'^) ^^"^ 
its coordination number z, defined as the average num- 
ber of interactions (contacts and distant attractions) for 
a particle in the packing, when the applied pressure is 
significantly smaller than Fq/o {P* ^ 1), as in the case 
of small powder samples assembled under gravity. With 
the values indicated above (at the end of Section lll Cp for 
toner particles, Fq/o? (the relevant pressure scale in 3D) 
is of the order of 100 Pa, which corresponds to a normal 
consolidation stress in a cohesive powder with 34% solid 
fraction [soj . 

In the absence of cohesion, the value of the applied 
pressure does not affect the properties of the packing 
(apart from setting the scale of intergranular forces) pro- 
vided the typical contact deflection, aPjKj^^ is small 
enough (rigid particle limit). We set this ratio to the 
value of F^jiaKti) in the cohesive case, i.e., equal to 
kT^ (see tablelU, so that typical contact forces are of the 
same order of magnitude (due either to P or predomi- 
nantly to Fq) in both cases. 

Effects of the initial solid fraction in the /-state, and of 
cohesion, friction and rolling resistance parameters on $ 
and z were measured in 3 sets of samples, with $i = 0.13, 
0.36 and 0.52. Each set consisted of configurations with 
the initial disorder (particle radii and initial positions) 
abiding by the same probability distribution, and the 
same number of particles (TV — 1400). The values of 
the friction coefficient, /i, used in these tests were 0.15 
and 0.5. The values of $ and z in these samples are 
listed in tables [llj and IIIII Each one is an average on the 
different samples, and the indicated uncertainly is equal 
to the standard deviation. 

Tables [ll] and IIIII show that the introduction of cohe- 
sion reduces the solid fraction at equilibrium, but this is 
a limited effect (less than 10% density reduction), which 
is quite insufficient to account for experimental obser- 
vations. Unlike powders or clays, 2D particle packings 
with $ > 0.7 cannot undergo very large plastic density 
increases. 

Theses tables also show that the increase of the friction 
coefficient and/or the inclusion of rolling resistance in the 
model tend to hinder motions and stabilize looser, less 
coordinated configurations, which results in a decrease 
of $ and z. 

However, the observed differences are rather small, es- 
pecially in cohesionless systems. To evaluate the infiu- 
ence of RR with a given value of /i, we define AX*^^) — 

(l-X^'^V^l^L): as the relative average decrease of the 
quantity X due to the existence of RR. For example, for 
/i = 0.15 results differ by a mere A^^'^ is) = 0.24% and 
^^(0.15) _ 1.3%^ and for /j, = 0.5, these variations are 
A$(0-5) = 0.84% and Az^-^) = 4.6%. Comparing the 
effect of RR on $ with /i — 0.5 for /i = 0.15 and /i = 0.5, 
one has A$(°-^V^*J''"°'^^'' = 3.5. Likewise, for coordi- 
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FIG. 2: (Color online) Aspect of force-carrying structures in cohesionless and cohesive samples. Contact forces are displayed 
with the usual convention that the width of the lines joining the centers of interacting pairs of disks is proportional to the normal 
force, on scale aP (left) and -Fo (right). Red, green, and blue lines distinguish compressive, tensile, and distant interactions in 
the cohesive case, while rattlers appear in grey in the cohesionless sample. 





Non-Cohesive samples 


no RR 


RR 


$1 


= 0.15 


M = 0.5 


= 0.15 


= 0.5 


Solid fraction 


0.130 ±0.001 


0.8262 ± 0.0007 


0.811 ±0.001 


0.8238 ±0.0014 


0.803 ± 0.002 


0.3631 ± 0.0006 


0.8256 ± 0.0005 


0.811 ±0.001 


0.8231 ±0.0013 


0.805 ± 0.002 


0.5244 ±0.0012 


0.8236 ± 0.0007 


0.8092 ± 0.0005 


0.8215 ±0.0005 


0.803 ±0.011 


Coordination number 


0.130 ±0.001 


3.174 ±0.012 


2.607 ±0.022 


3.160 ±0.024 


2.526 ±0.021 


0.3631 ± 0.0006 


3.187 ±0.025 


2.65 ±0.02 


3.123 ±0.013 


2.475 ± 0.025 


0.5244 ±0.0012 


3.181 ±0.015 


2.63 ±0.02 


3.15 ±0.03 


2.52 ±0.02 



TABLE II: Solid fractions and coordination numbers obtained at the preparation of the specimens in equilibrated samples 
under P/Kn = 10"^ for non-cohesive particles, using method 1. 



nation numbers z, one observes Az^°'^^ / Az^^-^^^ ~ 3.53. 
This shows a clear correlation of variations introduced by 
friction and RR. The data in the non-cohesive case also 
exhibit very little dependence on initial density $1. 

Results on cohesive systems show similar variations 
with the parameters of the contact model (friction and 
RR), but depend somewhat more sensitively on $1. 

More refined information on the contact network is 
provided by the distribution of local coordination num- 
bers, i.e. the proportions Xk of particles interacting with 
k neighbors, for < ^ < 6 (higher values were not ob- 
served). This distribution is depicted on Figure [31 for 
both cohesive and non-cohesive samples. These results 



gather information from all the statistically equivalent 
simulated samples, and slight corrections were applied 
in order to ignore the contacts with "rattler" particles 
in the non-cohesive case. Such particles are those that 
are free to move within the cage of their near neighbors, 
and transmit no force once the system is equilibrated. If 
they happen to be in contact with the backbone (i.e., the 
force-carrying structure) , then the forces carried by such 
contacts should be below the tolerance set on the equilib- 
rium requirement, and can safely be ignored. This is how 
the population of rattlers is identified. We observe that 
it can involve up to 18% of the total number of grains in 



the absence of cohesion (see Fig 2(a) I 
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Cohesive samples 


no RR 


RR 


$1 


= 0.15 


M = 0.5 


= 0.15 


= 0.5 


Solid fraction 


0.130 ±0.001 


0.7635 ± 0.0023 


0.751 ±0.001 


0.757 ± 0.002 


0.709 ± 0.001 


0.3631 ± 0.0006 


0.727 ± 0.001 


0.7232 ±0.0012 


0.710 ±0.002 


0.688 ±0.001 


0.5244 ±0.0012 


0.737 ± 0.002 


0.733 ± 0.002 


0.7248 ± 0.0002 


0.733 ± 0.002 


Coordination number 


0.130 ±0.001 


3.563 ± 0.005 


3.163 ±0.004 


3.189 ±0.008 


3.059 ± 0.003 


0.3631 ± 0.0006 


3.345 ± 0.009 


3.103 ±0.006 


3.253 ± 0.003 


2.971 ± 0.006 


0.5244 ±0.0012 


3.189 ±0.008 


3.059 ± 0.003 


3.096 ± 0.002 


2.851 ±0.001 



TABLE III: Solid fractions and coordination numbers obtained at the preparation of the specimens in equilibrated samples 
under P* = 0.01, Fo/{KNa) = 10^^ for cohesive particles, using method 1. 




FIG. 3: Distribution of local coordination numbers (percentage of total particle number), without (left graph) and with (right 
graph) cohesion. 



This contrasts with the cohesive case, for which nearly 
all the grains are captured by the force-carrying structure 
because of attractive forces, and the rattlers are virtually 
absent. The particles with one contact equilibrate when 
the deflection of that contact is ho, as defined in pT|) . 
With RR, such a particle is entirely fixed. Without RR, 
it is only free to roll without sliding on the perimeter 
of its interacting partner, because such a contact is able 
to transmit a tangential force smaller than or equal to 
IiKnIiq = hFq. 

Without cohesion, the coordination of the force- 
carrying structure can be characterized with a coordi- 
nation number z* , different from z: 



z* is the average number of contacts bearing non- 
negligible forces per particle on the backbone. Without 
cohesion, the backbone (or set of non-rattler grains) is 
the rigid part of the packing. With cohesion and RR, the 
whole interacting contact network is to be considered in 
order to study the rigidity properties of the system, and 
there are nearly no particles to eliminate. With cohesion 



and no RR, we observe in the samples obtained by the 
presently employed procedure (method 1) that the net- 
work of interparticle contacts or interactions is also rigid, 
apart from the free rolling of isolated grains with only one 
contact. (The rigidity properties of equilibrated samples 
arc discussed below in Section IIVI and Appendix |^ . 

Cohesive samples in equilibrium also comprise a small 
number of pairs of particles interacting without contact, 
i.e. separated by a gap smaller than the range of attrac- 
tion, Do. These are only a small fraction, below 1%, of 
interacting pairs. Such pairs do not contribute to dissipa- 
tion, since the frictional and viscous force components are 
only present in true contacts between neighboring grains. 
We observed that the time necessary to equilibrate the 
sample tend to increase when such distant interacting 
pairs are more numerous. 

In addition to the elimination of free rattlers, the most 
notable effect of cohesion on local coordination numbers 
(Fig. [3]) is to increase the proportion of disks with 2 con- 
tacts. Without cohesion, the Coulomb condition restricts 
the angle between the directions of the 2 contacts to val- 
ues between tt — 2<p and tt, where ip is the friction angle 
(tan</3 = ^). Thus if is small, a disk with two contacts 
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should have its center close to the line of centers of its two 
partners. The increase of the population of 2-coordinated 
disks as fj, is raised from 0.15 to 0.5 (see Fig. [3]) in co- 
hesionless systems corresponds to a less severe geometric 
restriction on contact angles. With cohesion, contacts 
may transmit a tangential force reaching iiFq while the 
normal force is equal to zero. Consequently, a disk might 
be in equilibrium with two contact points in arbitrary 
positions on its perimeter. As there is no geometric con- 
straint on the angle between the two contact directions, 
2-coordinated disks are easier to stabilize, and their pro- 
portion raises from about 5% without cohesion to above 
15% with cohesion in the case ^ = 0.15. A population 
of disks with one contact (therefore carrying a vanishing 
normal force, with deflection —h — ha) is also present. 
Those particles are fixed by a small rolling resistance, 
but are free to roll on their interacting neighbor without 
RR. Such a rolling motion is not damped in our model. 
Therefore, on waiting long enough, they should eventu- 
ally stop after a collision, in a stable position with 2 con- 
tacts. Such a collision is bound to happen because the 
contact network is completely connected. However, we 
stop our calculations when the kinetic energy is below a 
set tolerance (see Section III D[) . and we do not wait until 
all freely rolling disks reach their final position. Hence 
the remaining population of disks with one contact in 
samples without RR. 

The final configuration, with this preparation method, 
depends somewhat on the rate of compaction in the as- 
sembling stage. The latter is related to the choice of 
the dynamical parameter M, the "mass" with which the 
changes in cell dimensions are computed with Eqn. ([3]). 
The slight influence of the initial solid fraction, $1, also 
relates to such dynamical effects: a lower value of $1 en- 
tails larger colliding velocities, which favors larger final 
solid fractions. 

Although some of the aspects of the model (in par- 
ticular the homogeneous shrinking imposed through the 
periodic cell dimensions in a dynamical regime) do not 
correspond to experimental conditions, configurations of 
type 1 should be regarded as typical results of fast as- 
sembling processes, in which the particles are requested 
to balance the external pressure before stable loose struc- 
tures can be built. When the toner particles mentioned 
at the end of Section III CI are first fluidized, and then 
settle under their own weight, a rough estimate of the 
settling time, assuming particles are settling individually 
in air, and fall over distances of order 1 cm, is ~ 1 s. 
Fig. m with the value Tq ^ 10~^ s corresponding to such 
particles, shows that the duration of the "method 1" com- 
pression process is a few milliseconds. In practice, due 
to the presence of the surrounding fluid, the packing of 
a powder in a loose state by settling and compaction of 
an initially fluidized state is therefore considerably slower 
than this numerical process. 

In the next section, we consequently turn to the oppo- 
site limit, in which the external confining pressure is felt 
only after large, tenuous contact networks are formed. 



B. Method 2 

1. Numerical procedure. 

The second method to prepare numerical samples al- 
lows for aggregate formation before imposing an external 
pressure. Along with method 1, its different stages are 
schematically presented on Fig. H) 

The aggregation phenomenon plays an important role 
in the experimental preparation procedure of refs. [svllssj 
in which powder particles in a fluidized bed collide and 
stick to each other. Then they settle under their weight 
when the upwards air flow is abruptly shut off. The nu- 
merical method was designed to reproduce, in some ide- 
alized way, the final state of a set of colliding particles 
in the absence of external force fields. In the initial dis- 
ordered low-density configuration (the same /-state as 
in method 1), particles are now attributed random ve- 
locities drawn according to a Maxwell distribution, with 
mean quadratic velocity Vo- 

We performed systematic sets of simulations of disk 
packings with Vq = 9A8V* (see Eqn. Vq is thus 

large enough for the initial kinetic energy to overcome 
potential energy barriers in the process of aggregation. 
(The dependence of the final packing structure on this 
initial velocity of agitation, or "granular temperature", 
in the assembling stage in systems with small RR will be 
studied in Section HV B 6p . 

Once launched with such random velocities the parti- 
cles are left to interact and stick to one another within 
a cell of constant size, forming larger and larger aggre- 
gates, as appears on the image marked "B2" on Fig. [H 
Eventually, all particles are connected to one another by 
adhesive contacts, and reach an equilibrium position. At 
this stage, the two degrees of freedom of the cell are set 
free, and the stress-controlled calculation proceeds with 
cii — <J22 = (or P* = 0) until an equilibrium state is 
reached. This relaxation step does not lead to any re- 
arrangement of the contact structure, it only entails a 
very small increase of the solid fraction (hence the values 
slightly larger than <I>i given below). The final equilib- 
rium structure exhibits large density inhomogeneities, as 
apparent on F ig. [H which are characteristic of aggrega- 
tion processes [39| . and will be quantitatively studied in 
Section El 

Unlike cohesionless systems, which are devoid of any 
"natural" state of stress, clusters of cohesive particles 
can exist in a well defined state of mechanical equilib- 
rium in the absence of any external force. Once the state 
at zero pressure is obtained, we subsequently apply the 
same load P* = 0.01 as in method 1, which results in 
further compression and notable changes in the packing 
structure: $ increases from values close to $1 up to the 
0.45-0.55 range (see fig. E}. Nevertheless, the final solid 
fraction under P* = 0.01 is considerably lower than the 
one obtained with method 1. 

It should be noted on Fig. IH which summarizes the 
assembling procedures, that the aggregation stage makes 
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time (T^^ 



FIG. 4: Solid fraction versus time for both preparation procedures, showing some aspects of the configurations at different 
stages. Point A is the initial state (or $i). Aspects of configurations are shown for intermediate states Bi and B2, and for 
final equilibrated states Ci and D2 (at P* = 0.01). Point C2 corresponds to the stage when all disks are assembled in a unique 
aggregate, then equilibrated at P* = (both aggregation and equilibration stages take place between A and C2). The time 
unit is To = \J ma / Fo . Note the duration of the preparation process with method 2, and the difference in final equilibrated 
states compared to method 1. 



method 2 computationally quite costly because of the 
time necessary for clusters to merge, and especially for 
the stabilization of loose samples in equilibrium config- 
urations (lower contact numbers implying lower rates of 
energy loss as well as larger and slower fluctuations of 
soft, tenuous structures). In an attempt to limit the in- 
fluence of compaction dynamics, which results in denser 
samples when the lower density of the initial state al- 
lows the compaction process to accelerate more (as noted 
in Sec. IIII A[) . we tested the effect of limiting the max- 
imum strain rate Smax- Without any limitation, we ob- 
tained a maximum value e ~ 0.15 T^^. Using the sam- 
ples with <I>i — 0.13 (the lowest value used in this work) 
with Kn = lO^Fo/a, three different values for Cmax were 
tested: 0.10 T^^, 0.05 Tg"^ and 0.015 Tq-^ The condi- 
tion e > 0.10 Tq^ gave a final state close to the original 



one. The others two values produced similar results, with 
a relative decrease in density of about 10% compared to 
the original procedure. We chose to enforce condition 
Emax = 0.05 Tq"^' to saye computational time. This value 
has been applied to prepare all samples studied in the fol- 
lowing. 

Fig.[3]shows that method 2 succeeds in stabilizing open 
structures. Final solid fractions agree with the typical 
values observed in powders if one uses the correspon- 
dence between 2D and 3D packing fractions suggested 
by Campbell in [69j : 

$3D = 7^*?^' ^ 0.752$?/'. (14) 

Numerical samples under P* = 0.01, with solid fractions 
around 45%, would correspond to a powder consolidated 
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in the laboratory under 1 Pa with a sohd fraction of about 
23%. This is in satisfactory agreement with the experi- 
mental results of Ref . [5^ . 

We therefore regard method 2 as an appropriate way 
to reach an essential objective of this work, since stable 
loose structures are obtained. 

Although we perform simulations of a mechanical 
model, the final configurations exhibit at first sight 
(Fig. H]) similar features as those obtained with geomet- 
ric algorithms implemented in numerical studies of col- 
loid aggregation models [H, We are not aware of 
similar results in the literature, at least with equilibrium 
requirements comparable to those of Sec. IIIDI 

Tenuous, fractal-like contact networks contain denser 
regions and large cavities. Such heterogeneities produce 
long range density correlations, to be analyzed in Sec.lVl 
Without tensile contact forces, the walls of the cavities, 
comprising particles that are pushed towards the hole 
by the resultant of contact normal forces, would tend to 
buckle in. 

We regard method 2 as yielding typical results for as- 
sembling processes in which particles form tenuous aggre- 
gates before they are packed in a structure that is able to 
sustain a confining stress. In the sequel, we focus on the 
tenuous structures obtained with method 2. 



2. Global characterization of loose packings at P* = and 
P* = 0.01. 

We simulated four samples with 1400 disks and three 
of 5600 for $1 = 0.36, rather than lower initial densities, 
in order to achieve statistical significance at affordable 
computational costs, and to check for possible size effects. 
This set of samples will be denoted as series A. 

Samples with $ = 0.13 (series AO), which require the 
initial cell to shrink more before a stable network can 
resist the pressure, request longer calculations. Although 
some samples were prepared at P* = 0, we do not use 
them any more in the following, except for the values 
showed in Table |TVl 

To accelerate the numerical assembling procedure, we 
also created samples with = 10'*Fo/a, using an in- 
termediate value of $1 = 0.26, and softer contacts, such 
that K = 10^ in the initial aggregation stage (recall the 
time step is proportional to A/^/Kn)- Once equilibrium 
was reached with P* = 0, we slowly changed the stiff- 
ness parameter from k — 10^ to k — 10"*, and recorded 
the final equilibrated configuration. This procedure is 
about ten times as fast as the normal one, and generates 
similar structures and coordination numbers as series A 
prepared with the same /-state density. We shall refer to 
this set as series B. 

In table |TV] we list the corresponding results for solid 
fractions and coordination numbers. In such data we did 
not find a significant difference between the two differ- 
ent sample sizes, and therefore we did not distinguish 
between sizes in the presentation of statistical results. 



The tenuous networks obtained with method 2 collapse 
on changing the pressure: Table |V] gives the new values 
of $ and z after the compaction caused by the pressure 
increase from P* = to P* = 0.01. 

Structural changes between P* = and P* — 0.01 are 
shown on Fig. [51 which illustrates by means of four se- 
lected snapshots the mechanism of the closing of pores 
in a 1400 disks sample of series A. The first image cor- 
responds to equilibrium at P* — 0, and the fourth one 
to equilibrium under P* — 0.01. The two others show 
intermediate, out of equilibrium configurations during 
the collapse. One may appreciate how the denser re- 
gions grow and merge while pores shrink. Fig. [5] also 
makes it quite evident that the size of 1400 disk sam- 
ples is not very much larger than the scale ^ of density 
heterogeneities (typical diameter of large pores or dense 
regions, which will be studied in Sec. |V|. These sys- 
tems will exhibit large fluctuations in their mechanical 
properties: the rectangular shape of the final configura- 
tion displayed on Fig. [5] shows that the disorder is large 
enough for the mechanical response of the system to be- 
come anisotropic. Isotropy should be recovered in the 
limit of large sample sizes, L ^ 

Finally, Fig. [6] displays the histogram of local coordi- 
nation numbers (percentage of particles interacting with 
k others, < fc < 6), for the same samples as those of 
Tables HV] and El {fi = 0.5, $i = 0.36). It is remark- 
able that this distribution, in spite of the large differ- 
ence in sample geometries, remains rather close to the 
one observed in the denser packings made with method 
1 (compare P* = 0.01, fi = 0.5 case), just like global co- 
ordination numbers take very similar values in samples 
prepared with both methods (see Tables IIIII and |V| , in 
spite of the very different solid fractions. 

An essential conclusion of the present study is there- 
fore, for one given material, the absence of a general rela- 
tion between the density of a cohesive packing and its co- 
ordination number, in spite of previous claims [26i] . Both 
quantities are determined, rather, by the conjunction of 
micromechanical laws and sample preparation history. 



3. Effects of micromechanical parameters 

Adhesion should enhance the role of sliding friction and 
rolling friction, because the limiting values for tangential 
contact forces and rolling moments are both proportional 
to the elastic repulsive part of the normal force, N'^ ( 
\T,j\ < fiNfj, |f^J| < firNfj). Consequently, contacts 
with the equilibrium value ho of the elastic deflection for 
an isolated pair of grains transmit no normal force, but 
are able to sustain tangential force components as large 
as fj^Fo and rolling moments as large as firFo- Those 
values might turn out to be large in comparison to the 
typical level of intcrgranular forces under low external 
pressure (P* <^ 1). Therefore, even very low values of 
fi and /ir should affect the final structure of equilibrated 
packings considerably more than in the cohesionless case. 
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no RR 


RR 


$1 {z = 0) 


$ 


z 


$ 


z 


0.1301 ± 0.0003 (series AO) 


0.1303 ±0.0003 


3.197 ±0.002 


0.1304 ±0.0003 


2.656 ± 0.007 


0.2649 ± 0.0006 (series B) 


0.265 ± 0.001 


3.123 ±0.004 


0.2652 ±0.023 


2.963 ± 0.006 


0.361 ± 0.007 (series A) 


0.3616 ±0.0003 


3.1407 ±0.0016 


0.361 ±0.009 


2.660 ± 0.004 



TABLE IV: Values of $ and z on equilibrating configurations at P* = with /i = 0.5. Samples with $i = 0.26 correspond to 
ft = 10*, the others to k = 10^. 





no RR 


RR 


$1 {z = 0) 


$ 


z 


$ 


z 


0.2649 ± 0.0006 (series B) 


0.448 ± 0.006 


3.235 ± 0.003 


0.42 ±0.01 


3.085 ± 0.005 


0.361 ± 0.007 (series A) 


0.472 ± 0.008 


3.175 ±0.003 


0.524 ± 0.008 


2.973 ± 0.004 



TABLE V; Values of "I> and z in equilibrated configurations at P* = 0.01. These results are averaged over the whole set of 
samples prepared with \i — 0.5, for $i = 0.26 (series B) on the one hand, and for $1 — 0.36 (series A) on the other hand. Series 
AO, prepared with $1 = 0.13, yielded very similar results but due to computational costs the number of samples was too small 
to record data in statistical form. 




FIG. 5: Configuration of 1400 disk sample of series A without rolling resistance. Note the gradual closing of pores as the 
external pressure is increased from P* = (first image) to P* = 0.01 (last image) going through two intermediate stages. 
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Contacts per particle 

FIG. 6; Distribution of local coordination numbers in loose 
samples of series A obtained with method 2. Samples of series 
B gave a similar distribution. 



This is indeed the case for the coordination numbers ob- 
served in our simulations (see tables [IT] and IIIip which 
dropped more significantly, upon introducing the small 
level of RR we have been using, in cohesive systems than 



in cohesionless ones. 

On Fig. [7] we show the configurations at P* = (a) 
and P* = 0.01 (b) of the same sample assembled using 
method 2 with RR (left) and without RR (right). The 
denser regions in the inhomogeneous packings are joined 
by slender "arms" (see Figs. [S] and [7]). Such arms can in 
principle reduce to a chain of particles in the presence of 
rolling resistance. Such chains are otherwise destabilized 
by a rolling mechanism, hence the difference in the thick- 
ness of the arms with or without RR (see the blown-up 
detail in Fig.[7]-a), the lower coordination numbers of con- 
figurations assembled with RR (Tables HVl and IV)) . This 
might also explain the greater fragility of equilibrium 
configurations with RR, in which a larger compaction 
step (see Table IV)) is necessary, on applying P* = 0.01, 
before a new stable structure is reached. 

Another important parameter is the initial velocity of 
agitation, Vb- Its influence has been assessed on one 1400 
disks sample, with $1 = 0.36. The changes of coordina- 
tion number with Vb at P* = are presented on Fig. [51 

Low velocity values produce more tenuous aggregates 
[z ~ 2), since even a small level of RR is able to slow 
down local rearrangements and stabilize tree-like struc- 
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FIG. 7: Typical configurations of 1400 disk samples of series A with (left) and without (right) rolling resistance, at P* = (a) 
and P* = 0.01 (b). Note the difference in local structure of thin "beams" joining dense regions with or without RR. 



tures {i.e., devoid of flops) immediately after the colli- 
sions between particles or small clusters. 

A large kinetic energy cannot be absorbed by the RR, 
and as a result disks are able to rotate, which leads to 
better connected structures (z ~ 3). In a sense, a large 
Vq kills the effects of RR, and packings are similar to 
those made without RR in such cases. 

We therefore conclude that the connectivity of loose 




v/v* 



FIG. 8: Final coordination number z versus initial quadratic 
average velocity in agitation stage of method 2, normalized 
by characteristic velocity V*. The arrow points to the value 
most often used in our calculations. 



samples with RR assembled by aggregation depends on 
the initial magnitude of velocity fluctuations and on the 
level of rolling friction. 

As figure [5] shows, the same trend was found on reduc- 
ing contact stiffness parameter larger translational 
and rotational compliance creates more contacts. 

Vo is analogous to the particle fluctuating velocity in 
experiments o n g as-fluidized beds of xerographic toners 
under gravity [67[ . Such velocities are larger than the gas 
velocity by two orders of magnitude. Typically, one has 
Vgas ~ 1 — 4 mm/s, while V*, deduced from the contact 
parameters with relation ([7]) is about 1 cm/s. Such a 
value is therefore comparable to the particle fluctuation 
velocity. 

Of course, such a comparison is only indicative, be- 
cause the influence of Vq on packing structures depends 
on fir, and is also very likely to be affected to some extent 
by the viscous dissipation model we have adopted. Both 
rolling resistance and viscous forces are micromechanical 
features for which no accurate physical identification is 
available. Yet, it seems plausible that powder packings, 
because of their initial agitated states, stabilize in better 
connected states than predicted by geometric aggregation 
models. 

We now turn our attention in the next section to the 
forces in the contact networks, in particular the loose 
ones formed with method 2. 
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IV. MECHANICAL CHARACTERIZATION OF 
CONTACT NETWORKS 

Many numerical studies, in the past 15 years, have 
addressed the issue of contact network geom etry and 
force distribution in cohesionless systems [TOj. The im- 
age of force chains, i.e. a pattern in which larger inter- 
granular forces tend to line up on the scale of several 
grains, was evidenced in experiments [7l|, [T^I and sim- 
ulations [zl, [T^I, and the p.d.f of contact force values 
has often been measured and studied. An interpretation 
of the mechanical role of "force chains" [75] is that they 
carry the essential part of deviatoric stress, while the con- 
tacts carrying the lower forces are less sensitive to stress 
orientation and laterally stabilize the strong force chains 
against buckling. 

The main features of the distribution of forces and 
their spatial correlations have been reproduced by ap- 
proximate models (76j based on local equilibrium rules 
on each grain, supplemented by inequality constraints. 
One important such constraint is released in cohesive sys- 
tems, in which normal force components can have either 
sign. It is therefore worth investigating how the usual 
features of force-carrying structures in equilibrated gran- 
ular packings are affected by the presence of negative 
normal forces. One may also wonder to what extent the 
considerable difference in the density fields will affect the 
force patterns, given that the coordination of the force 
networks, as observed previously, does not seem to be 
very sensitive to density levels and density fluctuations. 



A. Force scale and force distribution 

The first obvious distinction between cohesive and non- 
cohesive systems is the appearance of a new force scale 
Fq, in addition to the one provided by the confining 
pressure, i.e., aP, the ratio of those characteristic forces 
defining the reduced pressure, P*. It is especially inter- 
esting to investigate the values and spatial organization 
of forces in systems with P* <C 1, as little information 
is to be found in the literature on this issue: numeri- 
cal studies of loose cohesive systems [2^ tend to focus on 
density and geometry of packings as a function of applied 
stresses. Some information on force networks is provided 
in a recent publication 11] on bead assemblies with cap- 
illary cohesion, but the confining stress is considerably 
higher is that study {P* of order 1) than in the present 
one. 

In the absence of cohesion, the distribution of force val- 
ues is usually presented in a form normalized by its av- 
erage, which itself scales with the applied pressure. This 
scaling can be made more quantitative on using a gen- 
eral relation between pressure P and the average normal 
contact force (Nij) and particle diameter d, which 

is known in the literature on powders as the Rumpf for- 
mula. We write it here in a form involving the spatial 



dimension D, which is valid both for D = 2 and D = 3: 

(15) 



„ 1 



In (|15p . d stands for the typical grain diameter. This rela- 
tion can be made more accurate if one notes that it stems 
from the standard formula for stresses in an equilibrium 
configuration (see the r.h.s. term in Eqn. [3|). To de- 
rive the formula, defining P ~ j) J2a=i '^aa, the average 
pressure, one assumes hij <^ Ri + Rj and then neglects 
correlations between particle radii and forces, assuming 



{N.,{R, + Rj))^FN{d). 



(16) 



Then, with a simple transformation of the sum, one ob- 
tains 



(17) 



With D = 2 and our diameter distribution (for which 
(rf) = 3a/4 and (d^) = 7a^/l2) this yields 



9 z$ 



(18) 



We found relation (|18p to be remarkably accurate in all 
our simulations, with or without cohesion, with configu- 
rations obtained by either method 1 or method 2, thereby 
checking that the correlations between particle sizes and 
contact forces could safely be neglected on writing (fT6|) . 

Without cohesion, Eqn. ([T8|) yields the correct scale for 
forces, i.e. the frequency of occurrence of intergranular 
forces larger than a few times F/v is very small. With co- 
hesion, when P* = or P* ^ 1, contact forces of order 
Pq are quite common, as shown on Fig. [51 on which nor- 
mal force distributions are represented. Hence Eqn. (jlSp 
cannot be used to predict "typical" contact forces. The 
presence of forces of order Pq explains the sensitivity of 
type 1 and type 2 samples with P* <C 1 to friction co- 
efficient and rolling resistance: densities and coordina- 
tion numbers (tables [1111 IIVI and W \ . in cohesive systems 
prepared under P* = or P* = 0.01 with /i = 0.15 
and with /i = 0.5, or with and without RR, differ sig- 
nificantly. Otherwise, if contact forces were of order of 
the average P/v , the value of which is correctly predicted 
by p^ . thresholds ^Pq or even ^r-Fb would be very large 
compared to typical forces and moments, and become ir- 
relevant. 

(It should be recalled that Rumpf 's name is often asso- 
ciated (as in Ref. [ll|) to a means to predict the macro- 
scopic tensile strength of a powder. As the essential in- 
gredient of the Rumpf approach [tJ is Eqn. [151 we refer 
here to that relation (like in [23]), as the Rumpf formula). 

Normal force distributions in cohesionless, cohesive 
type 1 and cohesive type 2 samples, the latter be- 
ing obtained with $i — 0.36 (series A), are shown on 
Fig. [21 Those distribution functions are roughly symmet- 
ric about 0, decay approximately exponentially at inter- 
mediate values, and vanish at — Pq, and Pq. In type 2 
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(c)With cohesion, method 2, P* = and P* = 0.01, fi = 0.5 

FIG. 9: Distribution of norm al forces for series A samples. 
The non-cohesive case (9(a) I is normalized by the average 



repulsive elastic part. The cohesive cases (9(b) and 9(c) \ are 
normalized by Fq (note that the average of the elastic part of 
is (Ne) ~ Fq in cohesive cases with P* <C 1). 



samples without RR, for P* = 0, there is a finite pro- 
portion of contacts carrying vanishing forces, about one 
fourth in the A series ($ = 0.36). In addition to this 
Dirac mass, there might be a power-law divergence near 
0, with an exponent our level of statistics is not sufficient 
to resolve accurately (about 0.6 to 0.8 in the range of 
forces between 10"'^ Fq and 10~^Fo). This proportion of 
zero forces is smaller, down to 9%, with RR, and drops as 
P* reaches 0.01, to 7% and 3%, respectively, without and 
with RR. It is worth pointing out that the correspond- 
ing contacts carry zero total forces, i.e. vanishing normal 
components {—h = ho, see ([TTt ) and no tangential elastic 
displacement either. In principle we cannot distinguish 
them from forces below the numerical tolerance defined 
in Sec. Ill Dl However, as we shall argue below in Sec- 
tion IIVBI under P* — one could expect all contact 
forces to vanish, and non-zero forces are related to the 
small, but finite degree of force indeterminacy. 

Before turning our attention to such features and to 
the spatial organization of forces, let us briefly discuss 
the differences between sets of (type 2) samples A and 
B. B samples, which are obtained with the "accelerated" 
procedure and $i = 0.26 (see Sec. IIIip . exhibit, due to 
their specific history, larger forces at P* = 0, with as 
many as 10% of the contacts transmitting normal forces 
N such that \N\ > Fq/IO, while this proportion lies below 
2% in A samples. On the other hand, B samples are 
looser, with more open contact networks under P = 
and a larger proportion of contacts (about one third in 
configurations without RR) carrying vanishing forces. In 
the following we shall use them to illustrate qualitative 
tendencies in very loose samples. 

When the pressure is increased to P* = 0.01, differ- 
ences in force distributions between A and B samples, 
despite their different sohd fractions (see table |V|, have 
considerably decreased, as shown on Fig. [TOl The influ- 
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FIG. 10: Comparisons between probability distribution func- 
tions of normal force values in samples of type A (histogram, 
in black) and B (shaded histogram, grey) without RR under 
P* = 0.01. 



ence of such differences in the aggregation stage as those 
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between our samples A and B are therefore expected to 
fade out after the systems are compressed to higher pres- 
sures and densities. 



B. Packing structure and force patterns 

The spatial organization of forces in type 2 samples, 
which we now discuss, is related to the distribution of 
force values, and should determine the ability of given 
configurations and contact networks to support stress in- 
crements. We first discuss systems without RR, then 
with the small RR values we adopted in most cases (see 
Table U]) . We emphasize the role of force indeterminacy 
and assembling history (the collisions by which cohesive 
clusters were built) in the final force patterns in equi- 
librium under vanishing or low applied stress. Extreme 
cases of systems with large RR on the one hand, or with- 
out friction on the other hand, are useful reference situ- 
ations, which we briefly examine and discuss. We con- 
clude this part with a discussion of the main physical 
implications of the relationships between force patterns, 
assembling process, geometry and micromechanical pa- 
rameters 



pattern of cohesionless systems, although there are of 
course compressive and tensile "chains". Unstressed re- 
gions are rather scarce in type 1 samples, although some 
areas with smaller forces are still present. The structure 
of type 2 samples under P* = 0.01 (Fig. [12]) is somewhat 
intermediate: isolated stressed clusters are still present, 
but elongated, force-chain-like structures emerge. 

To characterize such force patterns in a slightly more 
quantitative way, one can evaluate a threshold force Fperc 
such that contacts carrying a force F with ||F|| > i^perc 
percolate through the sample. Such a criterion was used 
to identify a "strong" subnetwork of force chains in [t^ . 
One observes Fporc of the order of the tolerance 10"'* Fq 
in $1 = 0.26 samples with no RR and P* — 0, which 
shows that stressed regions are isolated "islands" within 
the network. Fpoic raises to slightly less than O.IFq under 
P* = 0.01. 

Configurations of series A, assembled with $ = 0.36, 
possess the same qualitative features, although quantita- 
tively slightly weaker, due to their higher density. For 
instance, local stressed regions are somewhat less iso- 
lated, with a threshold force i^perc between 10"'^ Fq and 
8 X lO^^Fo at P* = 0. 



1. Qualitative aspects of force networks with no RR 

It is instructive to represent the forces carried by the 
contact network with a visualization of positive (repul- 
sive) and negative (attractive) normal forces, as was 
done on Fig. 2(b)[ showing the force network in one 
type 1 sample. Figs [11] and [T^ respectively correspond 
to equilibrated samples prepared with method 2 under 
P* = (immediately after the aggregation stage) and 
under P* = 0.01, without RR. They are represented here 
with (approximately) the same scale. Both belong to the 
(<I>i — 0.26) B series. Line widths, which are propor- 
tional to the intensity of the total interaction force, i.e. 
to ||F|| = 1 1 A^fi -I- rf| I , witness the presence, in spite of 
the low pressure, of many forces of order _fo (which cor- 
respond on the figures to line thicknesses comparable to 
particle radii). Stressed clusters, in loose type 2 samples 
under P* = 0, are separated by large parts of the inter- 
acting network in which contacts carry vanishing forces: 
the corresponding normal deflection is (Eqn. pT]) ) and 
there has been no elastic relative tangential displacement. 
Attractions (green) and compressions (red) have to com- 
pensate for Eqn. ([TS]) to hold true. This compensation 
appears to operate on a smaller scale in type 2 samples, 
because internal forces were previously balanced within 
isolated particle clusters. Such a local balance of forces is 
quite conspicuous at P* = (Fig. 111(1. in which internal 
stresses in small clusters often take the form of a periph- 
eral tension compensating a radial compression, or the 
other way round. This contrasts with samples prepared 
with method 1 (Fig. [2(b)]), in which the spatial distribu- 
tion of forces is more similar to the familiar "force chain" 



2. Force indeterminacy (without RR) 

The presence of large interaction forces of order in 
equilibrated samples is not obviously necessary, and is re- 
lated to the assembling process. Let us imagine particles 
are brought very slowly, one by one, within interaction 
range of the previous network, thus gradually building a 
unique cluster in equilibrium in the absence of external 
stress. One could expect, rather, each new contact to 
stabilize with N = T = and h = —hg. The existence 
of non-zero interaction forces in equilibrium is related to 
the hyperstaticity or force indeterminacy of the contact 
network. On writing all equilibrium equations for grains 
and collective degrees of freedom {i.e., setting accelera- 
tion terms to zero in Eqns. [T] [2]and[3]) and regarding all 
contact forces as unknowns, the degree of force indetermi- 
nacy H (or degree of hyperstaticity) is the number of re- 
maining independent unknowns, which cannot be deter- 
mined by the equilibrium requirement. If H = 0, knowing 
that some equilibrium forces exist (since an equilibrium 
state has been found) , then one would necessarily have all 
interaction forces equal to zero under P* = (since this 
is one obvious possible solution). The notion of force 
indeterminacy has been recently discussed by different 
groups in the context of granular materials, essentially 
because of the special case of rigid frictionless grains, for 
which the contact network is gencrically such that forces 
are uniquely determined [73,, 78, 79, .8Qj ■ The degree 
of force indeterminacy is linked to the number of degrees 
of freedom, equal to 3N (or 3A^ -I- 2 if the cell sizes can 
change), to the number of contacts Nc, the number of 
distant interactions A^^ and the number of independent 
mechanisms or floppy modes K (also called degree of hy- 
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FIG. 11: (Color online) Sample of type 2 (N=1400), in equilibrium under P* = after aggregation stage, with solid fraction 
$ = 0.26 (series B). Same conventions as on Fig. |2(b)| except for the blue color corresponding to contacts carrying a total force 
below tolerance 10~^Fo (deflection ho and no mobilization of tangential force). Note the large number of such interactions and 
the local compensation of attractions and repulsions in small prestressed clusters. To help visualize unstressed regions, disks 
only interacting at contacts bearing forces below tolerance are filled in light grey. 



postaticity [81[) by the following relation (written here 
for a fixed cell) 



3N + H^2Nc + Nd + K (no RR). 



(19) 



A proof of this simple result (which is classical in struc- 
tural engineering), and the relation of numbers H and K 
to the rigidity and stiffness matrices of the contact net- 
work, are recalled in Appendix [XI Mechanisms are those 
sets of velocities (or small displacements, dealt with as 
infinitesimal) which entail no relative velocities (or small 
relative displacements) in contacts. For distant interac- 
tions, only normal relative velocities are relevant, hence 
their particular treatment in (|19p . In Appendix |^ we 
explain how we determine whether a given configuration 



is rigid, i.e., devoid of mechanisms (apart from the two 
global translational motions of the whole set of grains, 
rendered possible by the periodic boundary conditions). 
It is customary to relate the level H of force indeterminacy 
to the coordination number z in granular materials. How- 
ever, this is not possible in general, which motivated our 
recalling (fTO|) in its complete form, can be rewritten, 
neglecting the very scarce distant interactions, as 

H = iV(.Z - 3) + K. 

Hence, in the absence of floppy modes, H = N{z — 3). 
However, there are still a few floppy modes on structures 
like those of Fig. [11] at the end of the aggregation stage, 
and this relation, which predicts a small degree of hyper- 
staticity relative to the number of grains (see Tables IIVI 
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FIG. 12: (Color online) Sample of Fig. 1111 with same scale and color conventions, in equilibrium under P* = 0.01. The solid 
fraction increased to $ = 0.39. The threshold force (used to distinguish blue lines and grey disks) was set to COIFq. 



and IV)) . is only approximate. Some mechanisms are due 
to the (exceptional) 1-coordinated disks and others, less 
trivial, are associated with larger parts of the structure 
which are connected to the rest of the packing via one 
single 2-coordinated disk. This floppiness is obviously 
related to the assembling process: before any external 
pressure is applied, nothing really requests the aggre- 
gates to possess a rigid backbone. The free motion of 
mechanisms in assembling method 2 is largely responsi- 
ble for the very long equilibration time (see Fig. 2]) : such 
motions entail no restoring force and no dissipation of 
kinetic energy. Floppy modes in the final state obtained 
with our criteria (Section III Dp being scarce (typically, 
a few such mechanisms per 1400-disk sample), we con- 
jecture that they would disappear entirely on adopting 
stricter equilibrium requirements in terms of kinetic en- 
ergy. If a mechanism survives, it should generically be in 
motion with a non- vanishing velocity, as a residual effect 
of the initially agitated state. As the connected aggre- 
gate partly folds onto itself, such motions should even- 
tually create new adhesive contacts, thereby reducing K, 
until the network becomes rigid. Once some rigid aggre- 
gate is formed in the assembling process, it will keep the 
same shape and structure, unless the collisions and per- 



turbations it subsequently undergoes cause it to break, 
because of the limited tensile strength of contacts or be- 
cause of the Coulomb inequality. This is the reason why 
the initial mean quadratic velocity of isolated grains in 
method 2 should be compared to F*, as given by ([7]). 

It is easy to see that the closing of one contact can con- 
vert an aggregate from floppy to hyperstatic, the simplest 
example thereof being the "double triangle" structure of 
Fig. [T31 By (fTO|) . this small structure, which is rigid 
(k = 3 counting the free motions of an isolated object in 
2D) has a degree of force indeterminacy H = 1. This is 
how the self-stressed clusters of Fig. [TT] are formed. Such 
structures have a strong influence on force values and 
force distribution. In particular, we show now that they 
entail specific correlations between normal and tangential 
force components in contacts. 



3. Local patterns and specific force orientations 

Fig. [14] shows all contact force values as points in the 
iV, T plane for a 5600 disk sample without RR equili- 
brated under P* — 0. Fig. [T3] displays a striking X- 
shaped distribution in the N, T plane, corresponding to 
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FIG. 13: Three (grey) disks initially forming an isostatic 
structure, when a fourth one (coming from the left) adheres 
to one of them (dotted position), it can roll (this is a floppy 
mode) until another (fifth) contact is formed, stabilizing it 
in the final position drawn with continuous line. The final 
"double triangle" structure is hyperstatic. Final forces (see 
main text) were computed for different initial velocities of the 
mobile disk and for different values of impact parameter S. 
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FIG. 14: Values of normal and tangential contact forces in 
a 5600 disk, type 2 sample, in equilibrium under P* — 0, 
with $ = 0.36 (A configuration). In addition to the remark- 
able cross-shaped pattern, marked with dashed lines of slopes 
±\/3, note the large number of very small forces, the numer- 
ous points with |T| <C \N\ and the relevance of the value of 
the friction coefficient (p = 0.5 here), as a small number of 
forces approach the Coulomb cone. 
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FIG. 15: Histograms of angle 6, between normal vector n and 
total contact force F. Conventionally, 6' = 0° for a repulsive 
normal force and T — 0, and 8 — 180° for a tensile normal 
force and T = 0. Shaded histograms (grey) correspond to B 
configurations ($i = 0.26), bold-line non-shaded ones (black) 
to A ones ($i = 0.36) 



a ratio T/N of ±\/3 This cross pattern fades away in 
systems which have rearranged to support P* — 0.01, 
although the corresponding specific T/N ratios are still 
overrepresented, as shown on Fig. [151 The cross pattern 
of Fig. [M] corresponds to angles 9 ~ 60° and 6 = 120° on 
Fig. [ini and the second graph shows that 9 ~ 120° still 
corresponds to a peak in the distribution once the sam- 
ple has been compressed (and rearranged) to P* = 0.01. 
As other characteristic features of force patterns in loose 
type 2 samples, this correlation between tangential and 
normal force components is stronger in the more tenuous 
networks of series B samples, for which the data are also 
represented on Fig. 1151 The difference between both sam- 
ple series tends to disappear on compressing to P* = 0.01 



(Fig.[T5(b)t. 

The prevalence of ratio | | — is in fact easy to un- 
derstand. Many disks are in equilibrium with two con- 
tact forces, with two other disks which are themselves 
contacting each other, as on Fig. (Occasionally, a 
third contact might be present, bearing a much smaller 
force, which we neglect in the present argument). In 
such a situation, without RR, the three equations ex- 
pressing the balance of forces and moments on disk D 
involve four unknown force components. Labels corre- 
sponding to contacts with disks marked 1 and 2 like on 
Fig. 1161 one obtains, on counting positively repulsive nor- 
mal forces on disk D and tangential forces with a positive 
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FIG. 16: The bottom disk, marked D, of radius R, is in con- 
tact with two other disks 1 and 2, themselves touching, whose 
radii are Ri and 7?2- At equilibrium, contact forces on disk 
D should be carried by the dotted line joining its 2 contact 
points, which determines the ratio of tangential to normal 
force components. 
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The ratio in ([^0]) varies for the radius distribution we are 
using in the present study but its most frequent value, 
corresponding to Ri = R2 = R, is VS. Fig. [17] shows 
the same graph as that of Fig. [HI in the case of a loose 
packing of disks with the same radii. In agreement with 
formula (|^D|) . the "X" shape is sharply defined. 




FIG. 17: Values of normal and tangential contact forces in a 
5600 disk, type 2 sample of monodisperse disks in equilibrium 
under P* = 0. Note the sharp "X" shape on blown-up detail 
of small forces. 



the contact corresponding to the common base of the 
two triangles should be very small, thereby explaining 
the "dense line" along the N axis on Fig. [TD It can be 
checked by direct inspection that local simple patterns 
like those of Figs. [T6| and [THj are indeed typical for the 
forces with ratios T/N around ±\/3, or with |r| ^ |A^|. 
The values of equilibrium forces within such a cluster 




FIG. 18: Hyperstatic 4-disk cluster, with 5 contacts. The 
force at the contact point between 1 and 2 should be carried 
by the continuous line joining this point to the intersection of 
the dotted lines. Those lines are respectively defined, as on 
Fig. 1161 by the two contact points Ci and C2 of lower disk D 
with disks 1 and 2, and the two contact points C'l and C2' of 
upper disk D' with disks 1 and 2. Note that the continuous 
line is close to the line of centers, hence a small value of ratio 
\T/N\ in the contact between 1 and 2. 

depend on how it was built. Without RR three disks 
forming a triangle equilibrate with zero contact forces, 
since there is no force indeterminacy. On simulating the 
collision of a fourth disk with such a triangle (as already 
sketched on Fig. [T3| all four particles having the same 
radius a/2, we could observe final equilibrium situations 
with contact forces depending on the impact velocity, 
provided of course a hyperstatic structure like that of 
Fig. [181 with 5 contacts, was assembled. Tensile forces 
equal to —0.133 x Fq in contacts Ci and C2 of Fig. [TH| 
were created for a contact with an initial velocity due to 
the sole acceleration of the distant attractive force over 
distance Dq = 10~'^a (within a range of impact param- 
eter S defined on Fig. [0)1 . Larger, repulsive forces were 
observed for higher initial approaching velocities. Self- 
balanced forces of order Fq therefore naturally appear in 
the assembling process. 



To understand the frequency of occurrence of very 
small T/\N\ values, let us now consider again the small- 
est cluster with force indeterminacy, without RR, which 
comprises four disks and 5 contacts, as schematized on 
Fig. [TSl Fig. [TSj shows graphically that the balance of 
contact forces implies that the tangential force within 



4- Systems with small RR 

With the small level of rolling resistance we have cho- 
sen, = 0.005a (see TableHl, the general features of the 
force patterns in systems without RR are only slightly al- 
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tered, as apparent on Fig. [191 which shows the interaction 
forces in an equihbrated sample of series B with RR un- 
der P* = 0. Like in type 2 systems devoid of RR under 
P* = 0, forces of order Fq only exist in isolated regions. 
Note, however, that the small forces outside these regions 
with self-balanced stresses do not vanish, but are of order 
HrFo/a, a feature which is further commented below. 

In principle, the discussion of force indeterminacy and 
rigidity is quite different with RR. Contacts now carrying 
a moment, the analog of relation \W\ becomes 

3N + H^3Nc + Nd + K { RR) (21) 

With RR (and as we could check with the method of Ap- 
pendix all connected clusters are rigid. One may 
therefore use directly k = (ignoring the two global 
translations) in (pij) . All independent loops contribute 
2 to the degree of force indeterminacy, and the coordina- 
tion munber corresponding to isostaticity (no Soppiness, 
no hyperstaticity) is equal to 2. To find self-balanced 
forces in a loop, note that each one of the two contacts of 
any particle in the loop will carry opposite forces (whence 
two independent force components); the resulting torques 
are then to be compensated by the rolling moments at 
the contacts, to be determined with a number of equa- 
tions (one per particle in the loop) equal to the number 
of unknowns (one per contact in the loop). However, 
those moments are severely limited by inequality (|10p . 
The constant force F transmitted around the loop should 
then be of the same order of magnitude as firFo/a , hence 
small. If we had used the same threshold for blue con- 
tacts on Fig.fTOlthan on Fig. [Til then all contacts within a 
loop, because they carry forces above the tolerance level 
IQ-'^Fo, would have appeared as red (compressive) or 
green (tensile). Resetting the threshold to 10~^Poj of 
the order of ^rFo/a, thus enabled us to distinguish the 
hyperstatic clusters analogous to the previous case with- 
out RR from the new source of hyperstaticity, the effects 
of which are limited by the smallness of the RR param- 
eter fir- We could check that the force threshold Ampere 
for percolation, as defined above in paragraph IIV B l] is 
close to O.OlPo in that case. If clusters made with RR, 
which are (infinitesimally) rigid could not be broken, no 
loop should appear because two independent clusters do 
not generically collide simultaneously in several points. 
The existence of loops in the final structure therefore 
witnesses the fragility of tenuous structures which form 
with a small level of RR (which is further confirmed by 
the large scale changes observed between P* = and 
P* = 0.01). 

Other features of force distributions and force pat- 
terns in systems without RR, such as the correlations be- 
tween normal and tangential contact force components, 
can still be observed with the small rolling friction level 
fir = S.lO^'^a. The graphs of Figs. [HI and [T51 if drawn 
for configurations prepared in the same way with a small 
RR, are very similar. The small RR level used in simu- 
lations therefore only introduces small quantitative dif- 
ferences in that respect, at least for the parameters of 



the assembling procedure defined in Section [III Bl In the 
next paragraphs, we investigate, first, as an instructive 
limiting case, the effects of large RR, and then the situ- 
ations in samples with low RR assembled with different 
initial random velocities (as on Fig. [8]). 

5. Effect of a large rolling resistance. 

Fig [20] shows the analog of Figs. [Tl] and [191 obtained 
with a large rolling resistance: /ir = 0.5a in a 5600 disk 
sample. The resulting structure has very few, large loops, 
hence an extremely small degree of hyperstaticity, and 
most contacts carry but very small forces. The charac- 
teristic prestressed clusters of Figs. [11] [T^l and [THl have 
disappeared. Such packings with large RR therefore ap- 
proach the limit in which a simple geometrical rule is 
adopted to aggregate particles: in the present case one 
recovers the results of the ballistic aggregation algorithm, 
stipulating that particles or clusters move on straight-line 
trajectories and join to form larger, rigid objects as soon 
as they touch. This results in isostatic, loop-free struc- 
tures with coordination number 2. The resulting contact 
network has no force indeterminacy, and is consequently 
not prestressed. Our introduction, in the previous sim- 
ulations, of a finite rolling resistance (and a finite fric- 
tion coefficient) changes those structures in two respects: 
first, they form better connected structures, with loops ; 
second, they carry significant self-balanced forces, of the 
order of the maximum tensile force in a contact. Those 
effects are however dependent on the initial conditions 
for aggregation, as we now report. 

6. Effects of initial velocities in aggregation process 

As shown on Fig. [51 the initial mean quadratic ve- 
locity Vq in the aggregation stage of assembling method 
2 determines the final coordination number of systems 
with RR. Isostatic, loop-free networks are formed with 
the small RR level {fir/o, = 0.005) used in our systematic 
simulation series provided Vb is small enough. The re- 
sulting force network, as displayed as an inset on Fig. [201 
approaches a tree-like, loop-free structure, in which all 
contact forces vanish under P = 0. 

Table IVII shows the dependence of coordination num- 
bers and force values on initial velocity parameter Vq/V* . 
One distinguishes three populations of contacts or inter- 
actions: those with repulsive, negative and vanishing nor- 
mal forces (i.e., below tolerance level 10~^Po)j and, like- 
wise, between the average number of contacts per grains 
of each kind, respectively contributing z+, Z- and zq to 
the coordination number z. A+ (respectively, A_) is 
the average value of repulsive (attractive) normal forces, 

N\_' (resp., N_') the quadratic average. These results 
illustrate the dependence of the force distribution on the 
initial velocity parameter. Force indeterminacy and sig- 
nificant non- vanishing forces appear as Vq/V* reaches 
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0.095 


2.004 


0.12 


0.12 


1.76 


0.046 


0.002 


0.047 


0.002 


0.95 


2.04 


0.38 


0.35 


1.3 


0.090 


0.002 


0.095 


0.002 


9.5 


2.66 


1.17 


1.23 


0.26 


1.7 


0.050 


1.6 


0.042 


95 


2.96 


1.46 


1.43 


0.07 


5.8 


0.16 


5.9 


0.096 



TABLE VI: Coordination numbers of repulsive, attractive and 
unstressed contacts, and values of the corresponding forces (in 
units of _Fo) in samples with RR prepared at different initial 
levels of agitation, as on Fig. [S] 



values of a few units, with VqIV* = 9.5 corresponding 
to the simulation series labeUed A and to the force dis- 
tribution shown on Fig. [Q] 

This set of results therefore bridges the gap between 
our mechanical studies of cohesive particle aggregation, 
with the parameters given in Tableland the preparation 



method of Section UlIBl involving parameter Vq, and the 
results of geometric algorithms, which are more tradi- 
tional in the field of colloid aggregation [40| . 

Geometric changes due to the breaking and rearrange- 
ments of clusters as they aggregate lead to better con- 
nected and presumably less fragile structures, which 
carry forces of the order of the maximum tensile force. 



7. The special case of frictionless disks 

As a complementary study of the opposite extreme 
case to that of large RR, we ran some exploratory sim- 
ulations of frictionless, cohesive grains (also devoid of 
RR). In the limit of rigid disks {k — *■ oo), one knows then 
that such assemblies are devoid of force indeterminacy: 
H = [zl, [tI, [13, HH . As a consequence, once large clus- 
ters are formed under no external pressure, all contacts 
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FIG. 20: (Color online) Same as Fig. 1191 in a sample with large RR, fir ~ 0.5a, N — 5600, and $ = 0.26. Inset: force network 
in A'^ = 1400 sample obtained with low initial mean quadratic velocity Vb and small RR (corresponding to the bottom left point 
on FigEI. 



should bear normal forces equal to zero. Such a situa- 
tion is depicted on Fig. [21] The aggregate represented on 
Fig. m] is obviously floppy. The analog of relation is 

2N + ii = Nc + Nd + K {fi = 0, no RR) (22) 

(It is customary, on counting degrees of freedom for fric- 
tionless disks or spheres, to discard rotations, which are 
all irrelevant, thereby reducing the number of degrees of 
freedom to 2N in the l.h.s. of ([22| ; an alternative is to 
regard each rotational degree of freedom as an indepen- 
dent mechanism). 

Formula (|22p in the frictionless case yields for H = 
a number of floppy modes equal to 2N — Nc — Nd = 
(4 — z)N/2. The configuration of flg. [^T] has a coordi- 
nation number z = 3.14, hence a number of mechanisms 



larger than 40% of the number of particles. Such ag- 
gregates are therefore very floppy, although particles are 
firmly tied to their contacting neighbors. Large parts of 
the particle cluster of Fig. [21] are connected to the rest of 
the structure by only one or two contacts, thereby allow- 
ing large scale motions maintaining all contacts. Not sur- 
prisingly, the application of a small pressure P* = 0.01 
to the system of Fig. [21] produces a very large compres- 
sion step, resulting in the configuration shown on Fig. 1221 
The coordination number is now 4.01 (corresponding to a 
very small degree of hyperstaticity, due to finite contact 
stiffnesses uKm/Fq = 10^, as well as to distant interac- 
tions), and the force-carrying network has a rigid subset. 

We conclude that assemblies of frictionless, cohesive 
particles are rather singular, and do not seem capable 
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FIG. 21: (Color online) Same as Fig. 1111 in a sample with 
iV = 1400, $ = 0.26, no RR and no friction (m = 0). 




FIG. 22: (Color online) Sample of Fig. [2T] under P* = 0.01. 
<1> increased to 0.72. 

of forming stable loose structures under a non-vanishing 
confining pressure. It could of course be conjectured, 
like in the frictional case, that floppy networks as shown 
on Fig. [5T1 with some residual motion, would gradually 
form better coordinated structures and eventually be- 
come rigid, but such an evolution is too slow to be ef- 
ficiently followed in our simulations. 



C. Discussion 

The study of force values, force distributions and spa- 
tial force patterns we have been presenting here opens 



quite a few perspectives that are worth pursuing in 
more detailed and quantitative form. In particular, we 
have left the investigation of elastic moduli and vibra- 
tional eigenmodes of the tenuous structures formed with 
method 2 for future work. 

However, two qualitative conclusions can be drawn, 
which might have broad physical relevance. 

First, essentially by direct inspection of force patterns, 
we observed that, in loose configurations under relatively 
low pressure if compared to the tensile strength of bonds 
(as expressed by P* <C 1), local arrangements of grains 
tend to form isolated self-stressed clusters where forces 
are of the order of the maximum tensile force in a con- 
tact, Fq. Those clusters comprise any number of grains 
between a few units to a few tens, keep the memory of 
the assembling process, and strongly influence the force 
distributions. These features are more apparent at lower 
densities. The degree of force indeterminacy H might be 
a useful indicator, but is not sufficient in itself, as it is re- 
lated to the coordination number, which is very similar in 
type 1 (dense) and type 2 (loose) systems, and, moreover, 
does not account for the role of inequalities ^ and p^ . 
As a general rule, loose cohesive systems tend to have a 
wider force distribution when H is larger, whereas the 
opposite behavior was observed for confined cohesion- 
less granular materials W^. Dense hyperstatic clusters 
in loose packings are connected by regions which bear 
very small forces. On increasing the applied pressure by 
small amounts, important changes occur, in which these 
prestressed regions merge together and large forces tend 
to organize in locally preferred directions, as in "force 
chains" . Such structures are likely to play an important 
role in the mechanics of loose cohesive granular assem- 
blies. Our mechanical study stresses the different effects 
of the two physical origins of forces - interparticle attrac- 
tion and applied pressure - which tend to create different 
geometries, force patterns and force distributions. 

Second, the structure of the loose packings and the 
forces they carry are strongly influenced by the assem- 
bling conditions. The relative duration of compression 
and aggregation processes might produce results as dif- 
ferent as type 1 or type 2 configurations. The velocity 
of agitation in the initial assembling stage affects the fi- 
nal coordination number, as shown in Sec. IIVB6I Such 
parameters affect the force patterns as well, and those 
are also modified if contacts are initially modelled as soft 
(k = 10^), as in the procedure leading to configurations 
B. 

We expect that mechanical strength properties will 
also be sensitive to the aggregation process. 



V. GEOMETRIC CHARACTERIZATION 

A. Introduction 

Aggregation processes are well-known to produce frac- 
tal structures, which have been studied for many years. 
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in particular with numerical simulations (see ref. [40| 
for a review). Universal fractal regimes due to vari- 
ous types of aggregation processes (ballistic, diffusion- 
limited, reaction-limited) are most conveniently observed 
in very low density samples. Indeed, an object of frac- 
tal dimension dp extending over distance L in c? dimen- 
sions (d > dp) will have an apparent volume fraction 
$ ~ L''^"'^, which vanishes as L ^ oo. Starting from N 
isolated particles in a finite volume with periodic bound- 
ary conditions, an aggregation process cannot produce 
a fractal geometry over arbitrarily large length scales. 
In practice, for low enough values of the aggregation 
process will begin just like in the $ — > limit, when cor- 
relations between particles can be neglected. Later on, 
the crowding and interpenetration of clusters will pre- 
vent the fractal behavior to extend to larger scales |83|, 
and a classical geometric model to describe this situation 
is a dense packing of fractal domains (sometimes called 
"blobs") of typical diameter ^. ^ is the upper limit of the 
fractal regime, and is related to <& (see the discussion of 
eqn. 1 in [40|) as 



(23) 



a relation which should be independent of the total sam- 
ple diameter L, provided L ^. This "fractal blob" 
model is reminiscent of semi-dilute polymer solutions [s^ 
and has been employed in many different physical situ- 
ations, e.g., silica aerogels [s^]- It has been shown to 
describe experimental results on the packing of cohesive 
powders [Ij, |5^ . If such a geometric description applies 
to our loose systems, then ^ should be of the order of 
the typical size of large density inhomogeneities (dense 
regions or holes) in the samples depicted on Fig. [Sj 



B. Definitions 

Self-similarity is conveniently detected on studying the 
density autocorrelation function (DACF), as follows. Let 
x(r) denote the indicator function of solid particles, tak- 
ing values 1 if point r is within a solid disk, and zero 
otherwise. Then we define the DACF as: 

C{r) = (x(R)x(R + r))R = ^ y x(R)x(R + r))dR, 

(24) 

with an average over the origin position R over the whole 
sample surface, of area A. On computing C(r) periodic 
boundary conditions should be accounted for, so that po- 
sition R + r stays within the simulation cell. Isotropy 
ensures that C{r) is only dependent on distance r — ||r|| 
in the large sample limit (or on taking its ensemble aver- 
age). C(r), by construction, takes the value (the solid 
fraction) for r = 0, and tends to as r ^ oo. 

In practice it is convenient to calculate, rather than 
C{r) — <I>^, its Fourier transform, a function of the mag- 
nitude k = ||k|| of wavevector k by isotropy, which we 



denote as I{k). I{k) is simply related to the Fourier 
transform x of the field x(r) by: 



I{k) 



A 



(25) 



The notation I{k) is of course reminiscent of the scatter- 
ing intensity per unit volume for wavevector k (as used 
in e.g. small-angle X-ray or neutron scattering experi- 
ments), which is equal to I{k), up to a "contrast factor", 
replaced by 1 in (|^ . 

A fractal structure with dimension dp va 2D should 
have a power-law decreasing scattering intensity over 
some range of k: 



^ a 



(26) 



An exponential cut-off of the decreasing power law be- 
havior of C{r) around ?■ '--^ ^ is sometimes used [85l. [86|: 



C{r) 



(27) 



where the length £, introduced to make C(r) appropri- 
ately dimensionless, is a constant of the order of the av- 
erage particle radius. Then the corresponding form of 
/(fc), is given in terms of Gauss's hypergeometric func- 
tion 2Fi[a,6;c;x] [83|: 



I{k) = Cst + $27r£2r(di. 



l+dp dp 



(28) 

In 2D, as soon as the particles form one continuous ag- 
gregate, the empty space is split into a set of disconnected 
holes or pores. The distribution of sizes and shapes of 
such holes is another way to characterize the system ge- 
ometry. 



C. Procedure 

To compute I{k) from the configurations obtained in 
simulations, we first discretized the density field x(r), 
i.e., we considered its values on the points of a regular 
mesh, with spacings Aa; and Ay along the edges of the 
rectangular cell of the order of a/100. x(k) was then 
evaluated using a two-dimensional FFT algorithm, from 
which I{k) was deduced by formula ([^S)) and orientation- 
ally averaged on binning values of wavevectors k accord- 
ing tofc = ||k||. 

The field x(r) — defines a set of holes. We charac- 
terize a hole labelled as H by the value of its equivalent 
radius Rh- Rh is defined as the radius of a disk with 
the same radius of gyration as the hole. Specifically, if 
Nh is the number of mesh nodes in the hole, which are 
labelled as i, 1 < i < Nh, and have coordinates Xi, yi, on 
denoting as (a;^,y^) the coordinates of the mass center 
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of the hole, one has: 



R 



2 



Nh 

E 



(29) 



Holes have complicated shapes and may be character- 
ized by other quantities such as eccentricity or higher 
geometrical moments, but such refinements lie outside 
the scope of this paper. 

For each sample, we record the first weighted moment 
(or mass average) of the distribution of hole equivalent 
radii, {R)w, defined as 



{R)w — 



J2h=i Rh 



(30) 



where n is the total number of holes in the sample. In 
loose cohesive samples, we obtained a rapid power-law 
decay for the shape of this distribution. Definition ((30)) . 
rather than a simple number average, ensures that the 
very small cavities (formed by three or four disks in con- 
tact) do not dominate in the evaluation of the average 
and {R)w indeed characterizes the large pores in the 
loose packings. However, this definition can only be ap- 
plied when holes do not percolate through the aggregate. 
Thus, we have restricted the calculation of {R)w to sam- 
ples with a non-vanishing confining pressure, P* > 0, in 
which case we regard it as an independent measurement 
of length scale 



D. Results 

Functions I{k) are shown on Fig. [23l along with their 
fits by Eqn. [Ml for P* = or 0.01 with and without 
RR, for the configurations of series A (parameters of Ta- 
ble [J and $1 = 0.36). The FFT calculations have been 
averaged over different characteristics density maps, and 
the bars denote the standard errors. To carry out these 
fits, we have applied the Levenberg-Marquardt method 
for nonlinear least-squares fittings [8^|. This fitting pro- 
cedure yields values oi dp and ^ listed in Table IVTIl As 
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P* = 0.01 


P* = 


P* = 0.01 
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1.925 ±0.024 


1.93 ±0.04 


1.53 ±0.04 


1.51 ±0.04 




8.29 ±0.15 


6.07 ±0.2 


9.3 ±0.4 


5.06 ±0.21 



TABLE VII: Fractal dimension and fractal blob size obtained 
on fitting the data of I{k) to Eqn. (pg)) . 

expected, the fractal dimension is conserved in the com- 
paction between P* — and 0.01, but the fractal range 
shrinks. The marked difference in dp caused by the intro- 
duction of a small level of RR is remarkable. While self- 
similar clusters are very nearly dense {dp approaching 
2) without RR, more open fractal structures are stabi- 
lized on small scales with = 0.005a. This value of the 
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FIG. 23: (Color online) Scattering functions I{k) of samples 
with and without RR for P* = and P* = 0.01, averaged over 
4 samples of 1400 disks and 2 with 5600 disks. Fits of data 
points with Eqn. [28] are drawn with continuous lines. Both 
with and without RR, I{k) is larger for P* = 0, corresponding 
to larger density fluctuations. 



fractal dimension obtained with RR appears to coincide, 
within the error bar, with the value dp = 1.55 ± 0.02 ob- 
tained for the ballistic cluster aggregation model [H, , 
assuming particles or clusters move on rectilinear trajec- 
tories and stick to one another, forming rigid objects, as 
soon as they touch. The tenuous, loop-free structure of 
such objects, as previously commented, is retrieved in our 
simulations on using a large rolling friction coefficient or 
a small level of initial velocity fluctuations. If measured 
on such samples as those of Fig. [201 the same result was 
obtained, as expected: dp — 1.56 ± 0.04. With small 
RR or larger initial velocities, our simulations produce 
structures with, apparently the same fractal dimension, 
but a larger coordination number. Another observation 
from Fig. [23] is the presence of a slight bump (maximum) 
in /(/c), for ^ ~ 10a, (which is not present in the fitted 
function [^5)1 . Such a feature is analogous to the peak 
in the structure factor of dense particle assemblies, and 
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is likely related to the packing of aggregates. As the 
aggregates are requested to be mechanically rigid they 
tend to be multiply connected and, at least in 2D, nearly 
impenetrable: the maximum in the structure factor is a 
signature of steric exclusion. 

The evaluation of the first weighted moment {R)w of 
the distribution of hole equivalent radii, for P* = 0.01 
yields {R)^/a = 6.6 ± 0.2, and (R)n,/a = 5.47 ± 0.14 
with RR. As expected, these results are similar to the 
values of ^ given in Table IVIII 



VI. CONCLUSIONS AND FINAL REMARKS 

Cohesive packings can form equilibrium structures at 
very small solid fractions in qualitative agreements with 
experiments in fine and ultrafine powders. Equilibrated 
configurations are sensitive to the level of applied pres- 
sure, relative to contact tensile strength Fq, as expressed 
by the dimensionless number P* . They crucially de- 
pend on the assembling procedure, as even a low pres- 
sure {P* <C 1) can lead to rather compact states if ap- 
plied to a initial "cold" (i.e., with vanishing or low ve- 
locities) granular gas of isolated particles, as in method 
1. If, on the other hand, particles are given some ran- 
dom motion and have time to stick to one another before 
having to sustain some stress, as in method 2, tenuous 
particle networks and open structures are obtained. The 
initial random motion, which is ballistic in our simula- 
tions, could be diffusive in practical situations in which 
fine particles are dispersed in a fluid. Some random rel- 
ative motion of different particles is also present, due to 
velocity fluctuations, in situations of global flow or sus- 
pension sedimentation. 

Under low pressure, such loose packings carry self- 
balanced forces of the order of the maximum tensile force 
-Fo in hyperstatic, well-connected lumps joined by thin- 
ner arms where many contacts carry vanishing or very 
small forces. Such structures are sensitive to the magni- 
tude of initial velocities with which particles collide on 
forming aggregates. In general force networks differ from 
the usual "force chain" patterns of cohesionless systems, 
and are associated with different force distributions. The 
force balance is strongly influenced by the structure of 
small aggregates that are first created on assembling the 
system. They evolve very fast as the system rearranges 
when P* grows even by small amounts (from P* = to 
P* = 0.01). 

Due to the limited strength of contacts with respect 
to tangential relative displacement and rolling, force- 
carrying structures therefore differ from the ones ob- 
tained with geometry-based algorithms in which any par- 
ticles or clusters that join form one unique rigid, unbreak- 
able object. The result of such algorithms is however 
retrieved, in the presence of rolling resistance, if large 
strength properties are attributed to contacts (to the RR 
parameter /i^ in particular) or if initial velocities of col- 
liding grains are kept low enough. In such limits isostatic. 



loop-free clusters are formed with coordination number 
2. 

Micromechanical parameters do otherwise influence 
the structure of packings and the initial (self-balanced) 
forces they carry, especially those without rolling resis- 
tance. 

The study of density correlations show that loose con- 
figurations can be regarded as dense packings of self- 
similar blobs of typical size ^ (about 10 times as large as 
the average diameter in our case), as in fractal clusters 
produced by colloid aggregation models. The estimated 
value of the fractal dimension, with RR, is compatible 
with the 2D result for ballistic aggregation, even when 
the connectivity (coordination number) is different. We 
thus expect different structures of the same density and 
fractal dimension to possess, due to the difference in loop 
numbers and self-stresses, different mechanical proper- 
ties. 

The fractal dimension appears to be larger in sys- 
tems without RR. Thus systems without RR seem to 
exhibit systematic qualitative peculiarities, and since a 
small level of rolling resistance is likely to exist in all re- 
alistic models, this feature should preferably be included 
in numerical studies. 

The effect of a growing pressure, as well as pressure 
cycles, on the packing density and internal state will be 
investigated in a forthcoming publication [42] . Other im- 
mediately related perspectives are the study of macro- 
scopic tensile and shear strength in relation to geometric 
characterizations and self-balanced forces. 



APPENDIX A: RIGIDITY AND STIFFNESS 
MATRICES 

Degrees of force indeterminacy H, of velocity indeter- 
minacy K and their relations are properties of the rigidity 
matrix G, which is defined as follows. First, let us denote 
as U a displacement vector for all degrees of freedom in 
the system, 

U = ((u,, A6'i)i<,<„, (ea)i<a<2) • (Al) 

in which one conveniently separates out in the displace- 
ment Ui of grain i the part due to the global strain, thus 
writing = — e • -|- Ui. U has dimension 3A^ -t- 2 for 
disks and 2 strain increments. Then for each one of the 
Nc contacts, say between i and j, the relative displace- 
ment of the contact point (with notations Rij for the 
radii, and tij for unit tangential vectors as in Sec. IIIB[) 

Suij = Ui + A9i X Riiij — Uj + A6j x Rjiij + e • , (A2) 

can be regarded as providing 2 coordinates to one 2Nc- 
dimensional vector of relative displacements, Su. As (jA2p 
expresses a linear dependence of Su on U, one has defined 
a 2Nc X (3A^ -I- 2) matrix, which is the rigidity matrix G: 

,5u = G • U (A3) 
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All coordinates of u and 6u are to be thought of as small 
(infinitesimal) increments, for which the system geome- 
try is fixed. The degree of displacement (or velocity) in- 
determinacy k is by definition the dimension of the null 
space of G. The relevant definition of relative displace- 
ments includes all relative motions that are associated to 
forces or moments. In the presence of RR, one should 
include all relative rotations S9i — 59 j into the compo- 
nents of 5u, the dimension of which thus raises to iNc- 
On the other hand, in the absence of friction the tangen- 
tial relative displacement of the contact point becomes 
irrelevant, and should only include normal relative 
displacements. In general all distant attractions between 
close neighbors should be dealt with similarly, because 
only normal forces are transmitted between such pairs. 
For future use we just denote as M the appropriate di- 
mension of the relative displacement vector. 

On writing 5u it is most convenient to use a local basis 
for each contact, with normal and tangential directions as 
coordinate axes. Increments of contact forces, and possi- 
bly moments (with RR), are related via the contact law 
to 5u. Together they define a contact force vector f , the 
dimension of which is equal to that of 5\i. f , in a system 
with RR, also includes rolling moments at contacts. 

Externally applied forces and torques onto the grains, 
as well as stresses, define together a vector of external 
forces F™': 

Y'^' = ((F„r,)i<.<Ar,(Aa„„)i<„<2). (A4) 

A denotes the surface area of the sample, so that the 
work of the load for small displacements is just F'^'^* • U. 
The equilibrium relations, stating that contact forces f 
balance the load F"'^*, just read (as one easily checks): 

poxt ^ . f ^ ^^5) 

with the transposed rigidity matrix, ^G. That matrices 
appearing in relations (|A3[) and (jA5[) are transposed to 
each other is just a statement of the theorem of virtual 
work: the work of external forces in any displacement 
vector is F'^''* • U = f • Su, provided F'^''* is related to f 
by (jASp and i5u is related to U by (|A3[) . By definition, 
the degree of force indeterminacy H is the dimension of 
the null space of ^G. 

The rank of matrix G is r = Nf — K, with Nf the num- 
ber of degrees of freedom (the dimension of displacement 
or external load vectors). This rank r is also the dimen- 
sion of the range of the matrix, which is the orthogonal 
subspace, within the M-dimensional space of relative dis- 
placements, to the null space of its transpose ^G in the 
dual space of contact forces. Hence r = M — H. We have 
obtained 

Nf + ii = M + K, 



which yields, according to the appropriate definition of 
relevant relative motions, relations (fTO]) . (|2T|) . and ([22| . 

Assuming elastic behavior in the contact {i.e., strict in- 
equalities in ([8]) and (fTO|) . which, as noted in Section llID| 
is the general case at equilibrium), in a quasistatic ex- 
periment contact force increments Af relate to relative 
displacement increments A^u with a contact stiffness ma- 
trix IC: 

Af = ^ • ASu. 

^ is a square, diagonal matrix, containing coefficients 
Kn, Kt and (with RR) Kr for each contact. Thus JC only 
contains positive elements, except for the (very scarce) 
distant interactions, which contribute the negative nor- 
mal stiffness —Fq/Dq in our model. If Af balances some 
load increment AF°'^', while 6u corresponds to the Nf- 
dimensional displacement vector U, one then has: 

^pcxt = g . U, 

where one has introduced the stijfness matrix K: 

K= "^E'^-E- (A6) 

(K is traditionally called dynamical matrix in the con- 
text of solid-statejDhysics and interactions of atoms or 
ions in a crystal :90] ) . Unlike G and ^G. K is always a 
square, symmetric matrix. It has to be positive definite 
in order for the equilibrium state to be stable, because 
it expresses the elastic energy associated with small dis- 
placements. (In fact, the full stiffness matrix also con- 
tains a small non-symmetric correction to (IA6[) [9l| due 
to the effect of contact forces prior to the application of 
the load increment, which we ignore here.) 

By construction, the null space of G is contained in the 
null space of K, and coincides with it in the absence of 
distant attractions, because K is then a positive matrix. 
In practice, the positiveness of K can be investigated with 
the Cholesky algorithm. We applied this method (in a 
form suitable for sparse matrices, stored in a "skyline" 
form) to the stiffness matrix of the contact networks of 
the simulated equilibrium configurations. This is how, on 
finding that K was positive definite, we could conclude 
that the contact structure was devoid of mechanisms (or 
floppy modes, eigenmodes of K with eigenvalue zero) in 
all cases with P* — 0.01. On the contrary, stiffness ma- 
trices associated with contact structures without RR at 
P* = usually possess some mechanisms, although we 
argued that their number k must be small. 
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